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HIGH-ORDER TWO-WAY ARTIFICIAL BOUNDARY CONDITIONS FOR NONLINEAR 
WAVE PROPAGATION WITH BACKSCATTERING* 

GADI FIBICHt AND SEMYON TSYNKOV*§ 

Abstract. When solving linear scattering problems, one typically first solves for the impinging wave 
in the absence of obstacles. Then, using the linear superposition principle, the original problem is reduced 
to one which involves only the scattered wave (which is driven by the values of the impinging field at the 
surface of the obstacles). When the original domain is unbounded, special artificial boundary conditions 
(ABCs) have to be set at the outer (artificial) boundary of the finite computational domain, in order to 
guarantee the reflectionless propagation of waves through this external artificial boundary. The situation 
becomes conceptually different when the propagation equation is nonlinear. In this case the impinging and 
scattered waves can no longer be separated, and the problem has to be solved in its entirety. In particular, the 
boundary on which the incoming field values are prescribed, should transmit the given incoming waves in one 
direction and simultaneously be transparent to all the outgoing waves that travel in the opposite direction. 
We call this type of boundary conditions two-way ABCs. In the paper, we construct the two-way ABCs 
for the nonlinear Helmholtz equation, which models a continuous-wave (CW) laser beam propagation in a 
medium with nonlinear index of refraction. In this case, the forward propagation of the beam is accompanied 
by backscattering, i.e. , generation of waves in the opposite direction to that of the incoming signal. Our two- 
way ABCs generate no reflection of the backscattered waves and at the same time impose the correct values 
of the incoming wave. The ABCs are obtained in the framework of a fourth-order accurate discretization to 
the Helmholtz operator inside the computational domain. The fourth-order convergence of our methodology 
is corroborated experimentally by solving linear model problems. We also present solutions in the nonlinear 
case using the two-way ABC which, unlike the traditional Dirichlet boundary condition approach, allows for 
direct calculation of the magnitude of backscattering. 

Key words, artificial boundary conditions (ABCs), two-way ABCs, radiation, the Helmholtz equation, 
nonlinearity, nonparaxiality, fourth-order schemes, self-focusing, backscattering 

Subject classification. Applied and Numerical Mathematics 
1. Introduction. The Helmholtz equation 

AE(xi, . . . ,xd) + k 2 E = 0 , A = + • • • + , (1.1) 
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models time-harmonic wave propagation in D dimensions. The simplest case is when k = ho, which corre- 
sponds to propagation of waves through a homogeneous medium. For example, in optics E is the electric 
field, k 0 = coono/c is the wavenumber, u> 0 is frequency, no is the (linear) index of refraction of the medium, 
and c is the speed of light. 

In many applications, one wants to solve equation (1.1) in the presence of an impinging wave field and 
boundaries, which can be either surfaces of obstacles or interfaces between different media. The impinging 
wave field is prescribed by a relation of the form 

E'mc = S; m pj n gj n g , (1*2) 

where Sj m pinging can, for example, be a point (specifies a spherical wave) or a plane (specifies a plane wave), 
and £P C is given. The physical properties of surfaces and/or interfaces, i.e., how they handle the impinging 
wave in terms of propagation through and/or reflection, are given by linear operator relations of the form 

L[E\ =0 on ^interface - (1.3) 

For example, if Interface is the surface of a perfect conductor, then (1.3) reduces to E = 0 on Interface (total 
reflection). 

Since equations (1. 1-1.3) are linear, one can solve the scattering problem in two sequential stages as 
follows. The solution is split into two components 

E — E\ nc + E s cat. • 

At the first stage one solves for the incoming wave field E inc , which is the solution of equation (1.1) in R D 
in the absence of any obstacles and/or interfaces, driven by the known source term (1.2). Typically, one can 
write this solution explicitly as a superposition of plane and/or spherical waves. Then, at the second stage, 
one solves for the scattered wave field E sc at , which satisfies equation (1.1) with no sources, subject to the 
boundary condition 


T[£scat] — T[£) nc ] On ^interface i 

which directly follows from (1.3). In the process of solving numerically for -E sca t> one has to replace 11° with 
a bounded computational domain. In doing so, one needs to introduce the ai'tificial boundary conditions 
(ABCs), see [28], which make the boundary transparent for outgoing waves and guarantee the solvability of 
the truncated problem on the finite computational domain, such that the computed solution is close to the 
original infinite-domain solution. 

In addition to the simplest case k = fc 0 , there are numerous applications where the medium is non- 
homogeneous, i.e., k = k(x ] , . . . , ;c /> ) . In this case, one may also need to solve for the incoming field E mc 
numerically (using ABCs), rather than analytically. However, as this problem is linear as well, one can still 
employ the linear superposition principle and thus first solve for E \ nc and then for .E scat . 

In the current study, we consider a more complex case when k depends also on the field intensity, i.e., 
k = k(wo, \E\ 2 ). For example, the propagation of an intense continuous-wave (CW) laser beam 1 through 
a Kerr-type medium such as water or silica, is described by equation (1.1) with k 2 = k q(1 + e\E\ 2 ), where 

*CW laser beam is a monochromatic wave, i.e., it is “purely” periodic in time, as opposed to, say, pulses and wave packets. 
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e = 4eocri2 and 712 in the Kerr coefficient 2 (e.g. [4,19]). In this case, beam propagation is governed by the 
nonlinear Helmholtz equation (NLH) 

AE+k 2 E = 0, k 2 = kl(l+e\E\ 2 ) . (1.4) 

Because of the nonlinearity, the equations for E- inc and E scat can no longer be decoupled as in the linear 
case. From a numerical point of view, this nonlinear coupling adds a new twist to the construction of the 
ABCs, since the Kerr medium interface at z = 0 is required to transmit E inc in one direction, and at the 
same time transmit E s cat in the opposite direction. Deriving and implementing this two-way ABC in the 
discrete nonlinear framework is a key emphases of this study. 

2. Physical model. Although our numerical approach is quite general, in order to motivate the pre- 
sentation we relate it to a specific physical problem, namely, that of an intense laser beam propagating 
through a nonlinear Kerr medium. The Kerr medium is located in the half-space z > 0, the directions of 
increasing and decreasing z are called right and left, respectively, and the wave source in the model is a 
right-traveling beam, impinging on the Kerr medium at z = 0. Therefore, the only physical boundary in 
the model is the transverse two-dimensional (x,y) plane at z = 0. For simplicity, we assume that the input 
beam is radially-symmetric in the transverse plane and denote the transverse coordinate by r = \J x 2 + y 2 . 

2.1. Two-way propagation of waves at media interface. At z = 0, the electric field E has both 
incoming and backscattered components. The value of the incoming wave upon entering the nonlinear 
medium is given by 

E mc (r , 0) = E? nc (r) . (2.1) 

In the current formulation of the problem, the two-way ABC at z = 0 has to ensure the reflectionless 
propagation of backscattered waves through the boundary (a radiation boundary condition ) and at the same 
time correctly prescribe the incoming signal (2.1). 

We note, however, that a more accurate physical model should include reflections from the media interface 
z = 0. These reflections can result in different values of the incoming wave field on two sides of the interface, 
i.e., E[ nc (r, — 0) ^ Ei nc (r,+0). In the current study we disregard this effect, which can be interpreted as 
either considering E® nc of (2.1) to be the part of the incoming wave that has already been transmitted past 
the z = 0 interface, or assuming continuity of the wavenumber across the interface. Similarly, we neglect 
the reflection of the scattered waves by the media interface at z = 0. In other words, we require that the 
boundary z = 0 be completely transparent for all left-propagating waves. In Section 8.2, we briefly comment 
on how one can incorporate a reflecting interface (i.e., discontinuity in k at z = 0) in the methodology that 
we are building. In fact, we consider this as one of the future extensions of our current work. 

2.2. Behavior as z — » +oo. Basically, as z — > +oo, we require that E have no left-propagating 
components. In this study we assume that at large distances propagation is diffraction-dominated and the 
field amplitude decays to zero, i.e., lim max \E(r, z)\ = 0, so that 

z — ^oo 0<r<oo 

lim k 2 = fcn . 

z — »+oo 

2 We note that the index of refraction is defined in the frequency domain. In the time-domain, the cubic nonlinearity 
becomes a nonlocal convolution, which , in the case of almost-monochromatic wavepackets, to leading order, is equal to a cubic 
nonlinearity [9]. 
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Therefore, at large z’s the solution is a linear superposition of right-traveling waves. 

In the discretization process (see Sections 5 and 6) we truncate the unbounded domain and introduce 
a far-held artificial boundary at z = z m ax . Similarly to the interface z = 0, the far-held boundary has to 
be transparent for all outgoing (i.e. right-propagating) waves. Consequently, the ABC at z = z max has to 
guarantee the rehectionless propagation of all waves traveling towards z = +oo. 

3. Paraxial approximation. Most research on wave propagation in a Kerr medium has been carried 
out in the framework of the nonlinear Schrodinger equation (NLS), rather than NLH. We now briehy describe 
how one derives NLS from NLH and quote some results on wave propagation in the NLS model. For more 
information on NLS theory, see, e.g., [11,19,26,27]. 

For reasons that would become clear later, we consider the NLH in R d with a general power-law non- 
linearity 

AE + k 2 E = 0 , k 2 = kl{l + e\E\ 2a ) . 


We denote the axial coordinate by z := Xd, and assume radial symmetry in the transverse plane of the hrst 
D — 1 coordinates, i.e. 

E = E(r, z) , r = \Jx\ + ■ ■ ■ + x 2 D _ x . 

We also separate the slowly-varying envelope ip from the fast oscillations and introduce nondimensional 
variables: 

V Z 

E = (r 0 koVe)~ 1/lT exp (ik 0 z)4’(r, z) , f = — , 5 = , 

ro zLdf 

where r 0 is the initial beam width and Ldf = 0o r o is the diffraction length. After dropping the tildes, the 
equation for the amplitude ip, in nondimensional form, is given by 

7n P ipzz + H’z + A±ip> + \ip\ 2,T ip = 0 , 


where the transverse Laplacian is 

8 2 d 2 _ d 2 D- 2d 

± dx\ dx 2 D _ x dr 2 r dr ’ 

and 

7np ( 2r 0 fc 0 

In typical physical setups the beam width ro is much larger than the wavelength A, which implies that 
0 < 7np 1 (or, equivalently, in dimensional variables, that ip zz <g; fco ip z ). Therefore, it is customary to 
employ the paraxial approximation , i.e., neglect the j n pip zz term. In that case, NLH reduces to the nonlinear 
Schrodinger equation (NLS): 

iip z + A ±ip + \ip\ 2a tp = 0 . (3.1a) 



The NLS is an evolution equation where z plays the role of “time” and the initial condition is given at z = 0: 


0(r,O) = £? nc (r) . 


(3.1b) 
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Therefore, under the paraxial approximation one approximates a boundary-value problem for the NLH 
with an initial-value problem for the NLS. Since the NLS accounts only for the forward-propagating wave, 
backscattering effects are neglected in this model. The question arises, therefore, whether and how the 
results of the NLS model remain valid at the NLH level, or alternatively, how these results are affected by 
backscattering. As of yet, almost no rigorous studies of these issues have been conducted. We therefore hope 
that the current study, which focuses primarily on developing a computational methodology for solving the 
NLH, will provide means for comparing numerically the NLH and NLS in the future. 

Let us now proceed with describing some specific results in the NLS model which are interesting to look 
at in the framework of the NLH. 

3.1. Critical self-focusing — arrest of collapse. We recall that the focusing NLS (3.1a) is called 
subcritical, critical or supercritical , when a (D — 1) is less than, equal to, or greater than 2, respectively. It is 
known that the solutions of both critical and supercritical NLS can actually develop singularities, i.e., blow 
up, at a finite z. There is, however, a marked difference between these two cases, as near the singularity 
nonlinearity dominates over diffraction in the supercritical case, while they are of the same magnitude in 
the critical case. As a result, unlike the supercritical case, singularity formation in the critical NLS is highly 
sensitive to perturbations, which can arrest the blowup even when they are small [11,12]. In this paper we 
focus on the critical case, which corresponds to the physical self-focusing (a = 1 and £> — 1=2). In that 
case, solutions of the NLS can become singular (i.e., blow up) after finite propagation distance, provided 
that their initial power ( L 2 norm) is above a certain threshold N c , which is called the critical power. 

The observation that the paraxial approximation breaks down near the singularity has been already 
noted by Kelley, in his celebrated paper on self-focusing [15]. Feit and Fleck [8] were the first to demonstrate 
that nonparaxiality of the beam can arrest the blowup, by showing numerically that initial conditions that 
lead to singularity formation in the NLS, result in focusing-defocusing oscillations in the NLH. In these 
simulations, however, they did not solve a true boundary-value problem for the NLH. Instead, they solved 
an initial-value problem for a “modified” NLH that describes the right-going wave only (while introducing 
several additional assumptions along the way). Akhmediev and collaborators [1,2] analyzed an initial- value 
problem for a different “modified” NLH; their numerical simulations also suggested that nonparaxiality 
arrests the singularity formation. Both numerical approaches ( [8] and [1,2]), however, did not fully account 
for the effect of backscattering. Fibich [10] applied asymptotic analysis to derive an ODE in z for self- 
focusing in the presence of small nonparaxiality. His analysis suggests that nonparaxiality indeed arrests the 
singularity formation, resulting instead in decaying focusing-defocusing oscillations. However, backscattering 
effects were neglected in this asymptotic analysis. 

Since there are no singularities in nature (i.e., the laser beam continues to propagate beyond the NLS 
blowup point), a natural question is whether initial conditions that lead to blowup in the NLS, correspond 
to global solutions of the corresponding NLH. To the best of our knowledge, the very issue of the solvability 
of NLH still remains unresolved, including the critical case a(D — 1) = 2. Therefore, we are interested in 
solving numerically the critical NLH as a true boundary-value problem , in order to address this question. 
Another issue of interest in the critical case is to calculate the amount of power which is backscattered for 
beams which do not blow up in the NLS model. We note that at present, there is no data coming from 
either analysis or numerical simulations, on the actual extent of backscattering, besides the general notion 
that it should be small. 

In order to simplify the calculations, we consider the critical NLH with £> = 2 and a = 2, i.e., 
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(3.2) 


^2+^\ E '(*’ r ) + # E = 0 > k2 = fc o(! + e l^| 4 ) > 

which corresponds to the critical NLS 

iip z + iprr + |V’|V = 0 - (3.3) 

Based on the insight gained from NLS theory, we can expect that the results for the critical NLH with D = 2 
and <7 = 2 would also apply for the critical NLH with D = 3 and a — 1. 

4. Nonlinear iteration approach. In this section we use a continuous formulation to outline and 
motivate the iterative numerical approach that we adopt in this study for solving the foregoing nonlinear 
wave propagation problem. The actual derivation, however, will be done completely at the discrete level in 
Sections 5 and 6. 

We are interested in solving the NLH (3.2) in the half-space z > 0, subject to boundary condition (2.1) 
for the incoming field, decay in the transverse direction 

lim E(r, z) = 0 , 

r— >oo 

and radiation conditions at z = 0 and z = +oo for the outgoing waves, as discussed in Sections 2.1 and 2.2. 
We build the iteration algorithm as follows. First, we define the linear version of the problem as 

L f [E} = 0, (4.1) 

where 

L F= fc2 + fr2 ^0 + eF ( r ' z )j > 

F(r, z) is a given function, and E satisfies the same boundary conditions as in the nonlinear problem 
we find the solution of the nonlinear problem (3.2) using the iterations 

L F(n) [E( n+ V] = 0, F n = \E ^\ 4 for n = 0, 1, 2, . . . ,Af , (4.3) 

with the initial guess E^°\r,z) = 0. Since there is no rigorous theory that guarantees the convergence 
of algorithm (4.3), our simulations (see Section 7) serve as a numerical test for the convergence of these 
iterations. In Section 8.3 we briefly discuss alternative approaches to the nonlinear iterations. 

4.1. Iterative solution of the variable-coefficient linear equation. In general, one can use any 

linear Helmholtz solver to solve equation (4.3) with respect to E^ n+1 ^ while keeping F^ frozen. In this 
study we solve (4.3) also iteratively as 

L 0 [£( ra+1 )] = -efcgirf") • £( m ) for m = 0, 1,2, . . . ,M(n) , (4.4) 


(4.2) 

Then, 


where 


r _ [ d 2 8 2 ] l2 

Lo ~ d^ + d^ +k °- 


Note that the function F^ does not change in the course of the iterations (4.4). 
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By rewriting formula (4.4) in the form 


E (m+1) = L -1 


- eklF <") • £?("*) , 


we see that it formally corresponds to the standard fixed point iteration scheme. Therefore, these iterations 
are more likely to converge when the RHS is small. We note that this occurs when eF^ <C 1, i.e., when the 
nonlinearity in the NLH is weak (fc 2 ss kf). We can expect this to be the case in physical self-focusing for 
the following reason. The Kerr coefficient of the medium n ,2 is so small that even for intense laser beams, 
upon entering the nonlinear medium, e|S? lc | 2 <C 1. In the framework of the NLS model, if the initial beam 
power is above the threshold for collapse, the nonlinear contribution to the index of refraction e|-E?| 2 (see 
(1.4)) would eventually become comparable to the linear one no- However, the asymptotic analysis in [10] 
suggests that nonparaxiality arrests self-focusing when €|I?| 2 < 1. As a result, k 2 ss kf for all z > 0. 

4.2. Direct solution of the constant-coefficient linear equation. At each iteration of the inner 
loop (4.4), we solve a linear constant-coefficient equation of the form 


L 0 E = $(r,z) , 


(4.5a) 


where the right-hand side (RHS) $ is given by 

$ = -ekfF^ • E (m) . (4.5b) 

Equation (4.5a), with $ given by (4.5b) and subject to the boundary conditions discussed earlier, is solved 
in the following way. We use Fourier decomposition in the transverse direction for the solution E, the RHS 
$, and boundary data Ef nc (r): 

E(r, z) = Y, ul (*) cos ( /r ) > $ ( r ’ -) = ^2 f( z ) cos (lr) , Ef nc (r) = ^ u° n ' c cos (Ir) . (4.6) 

i i i 

Because of the orthogonality of the Fourier modes, the l-th Fourier mode u l (z) of E(r, z) satisfies the ordinary 
differential equation 


u l zz {z) + k?u l (z) = f l (z) , kf = k 2 0 - l 2 , 


(4.7) 


subject to the Dirichlet condition for the right-going wave at z = 0 [cf. (2.1)]: 


«inc(°) = U i 


0,1 

inc 


(4.8) 


a radiation condition for the left-going wave at 3 = 0, and a radiation condition at z = +oo. It is at this 
level, i.e., after the separation of variables, that we implement the two-way ABC at z = 0 and the radiation 
boundary condition at 0 = + 00 . For that, we use the concept of the one-way Helmholtz equations . 3 

4.2.1. One-way Helmholtz equations and the radiation principle. Equation (4.7) admits two 
linearly-independent eigenfunctions: u W = and = e _ *v^ z . When kf > 0, u W = e*l fc, l z is the 

right-propagating wave and = e — I s i s the left-propagating wave, whereas when kf < 0, t/ 1 ) = 
is the right-decaying (evanescent) wave and is the left-decaying (evanescent) wave. Therefore, 

3 The term “one-way wave equation” is apparently due to Engquist and Halpern [7]. 
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the one-way Helmholtz equations that each admits only one of the two eigenfunctions while prohibiting the 
other one are: 


u z — = 0 , 

u z + = 0 . 


(4.9a) 

(4.9b) 


Equation (4.9a) corresponds to the right-traveling or right-evanescent wave u W, and equation (4.9b) corre- 
sponds to the left-traveling or left-evanescent wave u ^ . 

As mentioned in the end of Section 2.2, for the purpose of numerical solution we truncate the infinite 
domain [0, + 00 ) in 2 and reduce it to the finite interval [0, 2 max ]. The one-way Helmholtz equations (4.9) 
can be used as boundary conditions for equation (4.7) on the interval [0, 2 max ]. Indeed, if we want to make 
sure that near both edges of the interval [0, 2 max ] the solution is only composed of outgoing waves, then 
we need to use relation (4.9a) as the boundary condition at z = 2 max and relation (4.9b) as the boundary 
condition at z = 0: 


(4.10a) 

(4.10b) 


u z - iy kfu = 0 at 3 = 2 max , 

u z + iyjtfu = 0 at ^ = 0. 


Clearly, as the boundary conditions (4.10a) and (4.10b) each eliminate one of the two eigenfunctions 
and u^ 2 \ the homogeneous version of equation (4.7) on [0,2 max ] (i.e., when f l = 0) with these two boundary 
conditions is only satisfied by the trivial solution. Consequently, the non-homogeneous equation (4.7) with 
boundary conditions (4.10) is uniquely solvable for any RHS / concentrated on the interval [0, z max ]. From 
the standpoint of physics, the resulting solution is only composed of waves due to sources located inside 
[0, 2 max ], which radiate to the right and to the left, but contains no incoming waves from sources outside 
this interval. A solution of this type is said to satisfy the radiation principle. 

4.2.2. Adding the incoming power. As has been mentioned, for the particular problem that we are 
studying we also need to prescribe the incoming wave at 2 = 0, i.e., complement the radiation boundary 
condition (4.10b) for the left-traveling waves at 2 = 0 with a Dirichlet boundary condition (4.8) for the given 
right-traveling wave, which altogether will yield the two-way ABC. In the continuous framework, this can be 
done as follows. The incoming wave (4.8) gives rise to a solution of the form Substituting this 

expression into the one-way Helmholtz equation (4.9b), we arrive at the following inhomogeneous relation 


+ iy/^U = 2iy/^e^i^ e 


(4.11) 


As in the case of any inhomogeneous linear differential equation, the general solution to equation (4.11) can 
be written as a sum of the general solution v,u to the corresponding homogeneous equation (4.9b) and a 
particular solution u p to the actual non-homogeneous equation (4.11): 

u = «H + Up ■ 

We may pick the particular solution as the one generated by the incoming wave: u p = u^ c e z v^ z , and the 
general solution to (4.9b) is obviously given by «h = const • z . 

4.2.3. Obtaining the overall solution. In order to add the incoming power to the radiation solution, 
we replace the homogeneous boundary condition (4.10b) with relation (4.11) interpreted as a boundary 
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condition at the left edge of the interval: 


! kfu = 2 i 


K l “inc 


at z = 0 . 


(4.12) 


This implies that the overall solution will satisfy equation (4.7), subject to boundary condition (4.10a) at 
z = 2 m ax and boundary condition (4.12) at z = 0. Indeed, by linear superposition principle, the overall 
solution can be written as the radiation solution with the incoming power added: u = u ra diation + u^ c e tk,z , 
where u ra diation satisfies (4.7) and (4.10). A similar derivation in the finite-difference framework is presented 
in Section 6.5. 


4.3. Nested iterations. In summary, our solution algorithm consists of two nested iteration loops. 
On the outer loop (4.3) we perform iterations with respect to the nonlinearity for n = 0, 1, 2, . . . , J\f. On the 
inner loop (4.4) we solve the linear equation with variable coefficients (which we obtain at each nonlinear 
iteration) for m = 0, 1, 2, . . . , M{n). The numbers M = M(n) and A f, at which we terminate the inner and 
outer iteration loops, respectively, are determined experimentally in the course of iterations. 

Our particular choice of solver for the linear variable-coefficient equation (4.3) is motivated by the 
following two reasons: 

(I) The inner loop iterations (4.4) require inverting a linear constant- coefficient operator (which is 
the discrete analogue to L 0 ) rather than a variable-coefficient one. As a result, the inversion can be 
performed by a direct method that involves separation of variables and LU decomposition. Moreover, 
the implementation of the radiation boundary conditions, including the two-way ABC at z = 0, is 
particularly convenient to do with the operator £ 0 . 

(II) If we used some other linear Helmholtz solver, on each outer loop iteration (4.3) we would have had to 
invert a different linear operator Lp •*. However, using our particular linear solver involves a repeated 
inversion of the same operator throughout both inner and outer loops. This implies that the actual 
inversion can be performed only once in the very beginning and then the inverse operator, which is 
stored in memory, can be applied repeatedly to the changing right-hand side. From the standpoint of 
numerical efficacy this is beneficial because the inversion of the discretized L 0 amounts to performing 
LU decomposition of a family of sparse matrices obtained after the separation of variables. The result 
of the LU decomposition is also sparse, hence its application to a given right-hand side has only 
linear complexity. Since the number of iterations required for convergence is large (see Section 7), 
this yields substantial savings of computer resources. 

5. Discretization. We integrate the linear constant-coefficient equation (4.5) on a Cartesian grid of 
variables (r, z) in the finite rectangular computational domain [0, r max ] x [0, 2 max ]. Since the original physical 
domain stretches all the way to z = -l-oo, at the artificial boundary z = z max we set a radiation boundary 
condition that guarantees the reflectionless propagation of right-going waves (see Section 6). On the physical 
boundary z = 0 we set a two-way radiation boundary condition that similarly guarantees the reflectionless 
propagation of left-going backscattered waves and also correctly prescribes the right-going incoming signal 
(Section 6). As concerns the transverse direction r, we assume that the solution vanishes at r = r max : 


-L!(r ma x, 3 ) — 0 , z A 0 . 


(5.1) 


Physically, this condition amounts to having a conducting surface at r = r max , which acts as a perfect 
reflector. Therefore, we take r max sufficiently large so that reflections from this boundary do not contaminate 
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the solution in the primary region of interest near r = 0. We also assume that E is symmetric with respect 
to r = 0, i.e., 


E(r,z) = E(—r,z) , z > 0 . 


(5.2) 


This assumption is physically plausible and allows us to consider only half of the domain [0, r max ] in the r 
direction rather than the full domain [— r max ,r max ]. 

We use a uniform Cartesian grid with size h r and a total of M cells in the r direction ( h r = r max /M), 
and size h z and a total of N cells in the z direction (h z = z max / A' ) . Accordingly, the grid nodes are: 



r m =m ■ h r , z n — n * h z , m = 0, 1, 


M, n = 0, 1, 



(5.3) 


We discretize equation (4.5) using a fourth-order accurate central-difference scheme: 


KrE „ 


+ K Z Er, 


T /vq E rn n — 


m = 0 , 1 , . . . , M 


n = 2,3,. . . ,N 


(5.4) 


where 


L h rr E , n 


2,n T 16£/ m _i >n 30 E m n + 16 if 


m+l,n E m j r 2,n 


12 hi 


L z Z E m , ■ 


Em,n— 2 "b ^^E m n —i 30 E m n + 16£l min _)_i Em,n+ 2 

mf 


(5.5a) 


(5.5b) 


The index n that corresponds to the coordinate 2 runs from 2 to N - 2 in equation (5.4) because the 
stencil, which is five-node wide in each direction, obviously cannot be applied to any of the boundary nodes 
n = 0, 1, N — 1, and N located near z = 0 and 3 = 3 max . The treatment of these near-boundary grid nodes 
is discussed in Section 6 in the framework of the discrete radiation boundary conditions. 

Similarly, the direct application of the transverse part L^ r of the discrete operator in (5.4) may also 
require a special treatment of the near-boundary nodes m = 0,1, and M — 1. This treatment should take 
into account the transverse boundary conditions at r = 0 (5.2) and at r — r max (5.1). We can avoid this, 
however, by expanding the solution f? TO>n , for each n, in a finite series with respect to eigenfunctions of the 
transverse discrete operator L^ r , which also satisfy the two boundary conditions (5.1) and (5.2) [this is a 
discrete analog to the continuous Fourier expansion (4.6)]. This discrete eigenfunction expansion allows us 
to treat the operator L^ r in the transformed space from the very beginning and never implement it directly 
on the grid. In addition, the radiation boundary conditions in the 3 direction are most natural to implement 
in the transformed space separately for each longitudinal (i.e., ^-aligned) mode, as we have seen in the 
continuous formulation in Section 4.2.1. 

We shall now derive the discrete eigenfunction expansion for E m>n . Let us introduce the space of all grid 
functions that are equal to zero at m = M, i.e., 


V = 



m = 0, 1,...,M, ip M = 0 


}■ 


Clearly, for each n, the function E. n g V. We can define a weighted inner product on V: 


<t>) 


l 

2 M 


4’o(Po + 


1 M— 1 


m= 1 


(5.6) 
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Proposition 5.1. Let us consider a family of M one- dimensional grid functions of the argument m: 


'•Pm = cos ^(2fc - l)mA0^ , Ad = ^ ' k=l,2,...,M. (5.7) 

Then, 

(i) {f)}^cy. 

(II) The functions ipV) are orthogonal with respect to the inner product (5.6), i.e., 

(V> (fc ),V’ W ) = 0 fork^l. (5.8) 


(III) The set l x forms a basis in V. 

(IV) ^ are even functions of the argument m, i.e., symmetric with respect to m = 0: 

- 


(V) 


ipw are eigenfunctions of the transverse component of the finite- difference operator of (5-4) with 
eigenvalues A*, 4 , i.e., 


L h rr ^V = , 


A k — 


m 


16 sin 


2 1 ( - - sin 2 ((2 k - 1)A0) 


(5.9) 


if 


(k) 


m 


Proof. The inclusion (I) follows from the definition of the space V and the explicit form of the functions 
(5.7). To show the orthogonality (II), we calculate 


M - 1 


M — 1 


M ■ (ip( k \ipW) = ]T - \ = X) cos ((2 k - l)mA6) cos ((2/ - l)mA0) - * 

m= 0 m= 0 

i M_1 i 

= - ^2 [ cos ((2& + 2/ - 2)mA9) + cos ((2fc - 2l)mA6)\ - - 

m= 0 

i ^ i 

= - ^ [cos (2qmA6) + cos (2.smA0)] - - 

ra=0 
n M— 1 

= E' B 


p i2qmA0 _|_ ^—i2qmA0 _j_ ^i2smA0 _|_ i2sraA0 


ra=0 


1 _ e i2qMA9 ^ _ e ~i2qMA0 

+ 


J _ gi2qA6 _ g—i2qA0 


l 

+ 4 


1 — e 


%2sM AO 


+ 


1 — e 


-z2sMA0 ' 


^ _ gi2sA0 _ g— i2sA0 



We indeed obtain zero, because out of the two integer numbers q = k + l — 1 and s = k — l one is always odd 
and another one even, and thus one of the expressions in rectangular brackets on the last line in the previous 
chain of equalities is always equal to zero and another one is equal to two. Property (III) follows easily 
from the orthogonality (II) because the orthogonality implies that the M functions k = 1, . . . , M, are 
linearly independent, and the space V is obviously M-dimensional. Property (IV) is trivial and immediately 
follows from the definition (5.7). Finally, property (V), including the explicit expression for the eigenvalue 
A* given in (5.9), is obtained by directly applying the operator L f ) r of (5.5a) to each k = 1, ... , M. 


4 Note that for small wavenumbers the discrete eigenvalues and eigenfunctions are similar to those in the continuous formu- 
lation (cf. (4.6) and (4.7)) as A k as (k - l/2) 2 (7r /r max ) 2 and i/> \ T=mhr = 4>m ' = cos ((k - l/2)7rr/r ma x). 
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The application of L k r to a ip^ in the near-boundary nodes requires using the symmetry property (IV) and 
also noticing that all il^ k \ k = 1,... , M, are, in fact, antisymmetric with respect to m = M, which again 
immediately follows from the definition (5.7). ■ 

Proposition 5.1 shows that the system forms an orthogonal basis of the space V, composed of 

the eigenfunctions of the operator L\ f r , which are symmetric with respect to m = 0 and vanish at m = M. 
For all n we can construct the expansion with respect to these eigenfunctions according to 

1 1 M_1 

Uk,n = (E.,n,ip (k ^) = + Jf £ R m,n cos ((2fc - l)mA0) , k =1,2,... ,M , (5.10a) 

m= 1 

so that 


M 


M 


E n 




^k, 


cos ((2k — 1)toA0) = 2 £ Uk,n4' ( m ] , m = 0, 1, . . . ,M 


(5.10b) 


fc=i 


k=i 


Representation (5.10b) can be easily verified by directly substituting Uk, n of (5.10a) and performing the 
transformations similar to those performed when proving Proposition 5.1. Obviously, formulae (5.10a) and 
(5.10b) are particular versions of the direct and inverse discrete Fourier transforms, respectively. 

The above eigenfunction expansion can be used to implement the transverse discrete differentiation along 
with the boundary conditions at r = 0 and r = r max . Indeed, if we expand E m%n and the RHS d> m ,„ in the 
form (5.10b) with the coefficients Uk, n and respectively, obtained using (5.10a), then, because of the 
orthogonality of the eigenfunctions t(4 k ) (5.8), we arrive at the following family of one-dimensional discrete 
equations: 5 






k,n 


^k,n—2 H" 30 + 16?/fc, n +i Ufa' 


n+2 


12 hi 


+ k c Uk, n — fk,n > 


k 2 c = k%- A*, k = 1,2, ...,M , n = 2,3,... ,N — 2 , 


(5.11) 


where the eigenvalues {A*} are defined in (5.9). Each of the M equations of (5.11) is independent of the 
others and will be solved separately using the methodology of Section 6. Having obtained the modal solutions 
Uk, n for all k = 1, 2, ..., M, we then recover the overall solution E mn by means of the inverse transformation 
(5.10b). 

5.1. Implementation of transformations (5.10) using FFT. It is convenient to implement the 
direct and inverse transformations (5.10a) and (5.10b) using the standard discrete Fourier transform, for 
which library subroutines optimized for performance are available (fast Fourier transforms). To do that, we 
note again (see end of the proof of Proposition 5.1) that representation (5.10b) allows us to extend E m>n for 
any n beyond m = 0 and m = M using the explicit form of the basis functions see (5.7). The extension 
for negative to’ s is symmetric with respect to to = 0, and the extension beyond m = M is antisymmetric 
with respect to m = M. For a given function E m<n , m = 0, 1, . . . , M, it is convenient to extend it first anti- 
symmetrically with respect to to = M (so that the function be defined for m = 0, 1, . . . , 2 M), and then also 
extend it symmetrically with respect to to = 0 (so that it finally be defined for to = —2 M, ... , 0, . . . , 2 M). 

5 Note the analogy to (4.7). 
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In doing so, we arrive at a periodic grid function with the period AM. It is easy to see that for a function 
extended in this particular way the standard discrete Fourier transform 

2M — 1 

«i,„ = ^ E m , n e- ilmA0 , l = —2M , . . . , 2M — 1 , (5.12a) 

m=—2M 

reduces to (5.10a). Indeed, as E m ^ n is real we will always have = U-i, n , and in this particular case the 
symmetry with respect to m = 0 implies that all tq >n are also real and thus ui >n = U-\ %n . Consequently, we 
can consider only 2 M + 1 independent real coefficients ui n for l = 0, 1, . . . , 2 M. Then, the antisymmetry 
with respect to m = M will yield that u^ n = 0 for all even l = 0, 2, 4, . . . , 2 M and we are thus left with only 
the coefficients tq, n for odd l = 1, 3, 5, . . . , 2 M — 1. In other words, we can rewrite (5.12a) as follows 

1 1 M ~ 1 

= < 2 j[f Eo,n 1 /* ' E m n cos (IfiiA^) , l = 1,3,..., 2 M — 1 , 

m= 1 

and conclude that it indeed coincides with (5.10a) if we change notations from l = 1,3,5,... ,2 M — 1 to 
k = (l + l)/2, k = 1, 2, . . . , M. Similarly, it is easy to see that because of the aforementioned properties 
of ui^ n (ui tn = u/, n real, and u^ n = 0 for / = 0,2, 4, . . . ,2 M), the standard inverse discrete Fourier 

transform 

1 2M-1 

Em ’ n = JM ^ u Une ilmA \ m = — 2M, . . . , 2M , (5.12b) 

1 l=-2M 


reduces to (5.10b). 


6. The one-dimensional discrete Helmholtz equation. In this section we analyze the discrete one- 
dimensional linear hon-homogeneous Helmholtz equation (5.11), paying special attention to the treatment 
of the boundary conditions for z = 0 and 2 = 2 max . We recall that the boundary conditions at z = z max 
should guarantee that this boundary be transparent for all waves traveling to the right (i.e., a standard 
radiation ABC). The boundary conditions at z — 0 should guarantee that this boundary be transparent for 
all backscattered waves traveling to the left, and at the same time impose the given incoming wave field 
( two-way ABC). We emphasize that we have not discussed a particular discrete form of these boundary 
conditions until now, since typically the ABCs are most convenient to set in the transformed space rather 
than original space [28]. 

To simplify the notations, we drop the subscript k , so that equation (5.11) takes the form 


Un—2 T 30 Uji T ldUyj-i-i U n ^-2 

12/fr 


+ k 2 c u n 


fn, n 


2,3,... ,7V -2. 


(6.1) 


Equation (6.1) is a fourth-order difference equation. It is obtained, however, as a fourth-order accurate 
difference approximation to the second-order differential equation. Therefore, compared to its original con- 
tinuous counterpart, the difference equation (6.1) requires additional boundary conditions. A total of four 
boundary conditions are needed to guarantee the solvability and uniqueness for equation (6.1). Two extra 
boundary conditions that are not present in the continuous case are a pure numerical artifact. They are 
accounted for by the presence of two extra evanescent waves among the solutions of the homogeneous version 
of equation (6.1) in addition to the two standard traveling or evanescent waves (see Section 6.1). Altogether, 
these four boundary conditions should ensure the desired behavior of the solution near z = 0 and near 
- = -max- We also reiterate that the finite-difference equation itself obviously cannot be written in the form 
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(6.1) for the grid nodes n = 0, 1,1V — 1 and N. A special form of the discrete equation for these four grid 
nodes is therefore required; this special form will actually constitute the boundary conditions and make the 
total number of equations in the linear system be equal to the number of unknowns. 


6.1. The discrete homogeneous problem. We start by analyzing the homogeneous counterpart to 
the finite-difference equation (6.1) over an infinite grid domain, i.e., 


Un— 2 T ldu^—i 30 U n A 16u n+ i U n -\-2 


+ k 2 u n — 0 , 


n = 0, ±1, ±2, . . . 


(6.2) 


Proposition 6.1. Let a = ( h z k c ) 2 be such that either 0 < a < 16/3 or — 3 < a < 0. Then, the general 
solution to equation (6.2) has the form 


Un = Cl</( + C 2 «” + c —lQl " + C_2«2 " . (6-3) 

where ci,c 2 ,c_i, and c_ 2 are arbitrary constants, and qi and q -2 are roots of the characteristic equation that 
corresponds to (6.2). 

In addition, 

(I) When 0 < a < 16/3, q r ( and qf n are waves propagating to the right and to the left, respectively. In 
particular, when 0 < a 1, then 

qi = e ik ‘ h * + O ((k c ■ h z f) , (6.4a) 

q 2 = e ~ik.h z + 0 ( {kc . hz)6 j ? (6.4b) 

and as such, q " and qf n are the discrete analogues of the right and left traveling waves e lkcZ and 
e~' lkcZ , respectively, with fouHh-order accuracy. 

(II) When -3 < a < 0, q? and qf n are evanescent waves decaying to the right and to the left, respectively. 

(Ill) In both cases, i.e., for 0 < a < 16/3 and for -3 < a < 0, qtf and q% n are evanescent waves decaying 
to the right and to the left, respectively. 

Proof. Let us introduce the characteristic algebraic equation 

-1 + 16 q + (12a - 30)g 2 + 16g 3 - q 4 = 0 (6.5) 

for the homogeneous finite-difference equation (6.2). It is generally known (see, e.g., [14]) that if all the roots 
qj of a given characteristic algebraic equation are distinct, then the general solution to the corresponding 
homogeneous finite-difference equation is obtained as a linear span of the grid functions qf, where the power 
n is determined by the grid location. In the specific case that we are studying equation (6.5) is a quartic 
algebraic equation and thus provided that its four roots {q :/ }] =1 are distinct, the general solution to the 
homogeneous equation (6.2) has the form 

Un = Cl q r ( + C 2 ^2 + c 3l3 + C4<?4 » ( 6 -6) 

where {cj}j =1 are arbitrary constants. 

Hereafter, we restrict ourselves only to the case when the roots {</y }y =1 of equation (6.5) are distinct. 
By explicitly calculating {qj}j =1 (see below), we will show that multiple roots are only possible for the two 
cases a = 0 and a = 16/3, which are easy to avoid in practical computations. 

To simplify the actual calculation of the roots of quartic equation (6.5), we first note that by dividing 
(6.5) by q 4 we arrive at exactly the same equation for l/q. Therefore, if q is a root, then q~ 1 is also a root 
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(this follows, of course, from the fact that the discretization (6.1) is symmetric). Accordingly, we rename 
the four roots of equation (6.5): gi, < 72 , gj -1 , and </.J 1 , and write: 

- 1 + 16g + (12a - 30)g 2 + 16g 3 - g 4 = -(q - qi)(q - g^Xg - g 2 )(g - ^ 

— ( q 2 — dig + l)(g 2 — 6?2<7 + 1) = —1 + (d\ + d 2 )g — (2 + did 2 )g 2 + (di + c? 2 )o |3 — g 4 » 

where 


di — qi + q 1 1 , d 2 — g 2 + ?2 1 • (6-8) 

By comparing the beginning and the end in the chain of equalities (6.7) we obtain the following system of 
equations for di and d 2 : 


d\ + d^ — 16 , — 2 — did 2 — 12a — 30 , 


from which we find that 


d\ — 8 — 6-\/l + a/3 , d 2 — 8 A 6-\/l + a/3 . 


(6.9) 


From formulae (6.9) we conclude that both d\ and d-2 are real provided that a > —3. If, for example, 
h r « h z (the cell aspect ratio of the discretization is close to one), then the definition of k c (see (5.11)), 
where A * is given by (5.9), along with the definition of a = ( h z k z ) 2 , suggest that even for negative a’s their 
absolute values are sufficiently small and thus we can always assume that a > -3 and consequently, consider 
di and d 2 real. However, allowing for the complex values of d\ and d 2 may only make the analysis more 
cumbersome, but does not change any of the results hereafter. This, in particular, is corroborated by the 
computations of Section 7.1, which were conducted on the grids with cell aspect ratios 20/1 and 20/3. 

From (6.8) we have that qj and qj 1 are the roots of the quadratic equation 

q 2 — djq + 1 = 0, j = 1,2 . (6.10) 


Let us analyze the case j = 1 first. For 0 < a < 16/3, equation (6.10) has two complex conjugate roots 


Qi 


di + i — d 2 
2 



iy/ 4 - d\ 
2 


( 6 . 11 ) 


From (6.11) it follows that |<Zi| = 1^ 4 | = 1 and, in addition, that when 0<a€l then (6.4) holds. 
When — 3 < a < 0, we have 


Qi 


di~ „-i 

O ’ "l 


d\ + \J di\ — 4 
2 


(6.12) 


Therefore, both roots are real and satisfy |gi| < 1 and l^j -1 1 > 1, showing that g" and gf” of (6.12) are discrete 
analogues of two evanescent waves. We note that as a changes from positive to negative in formulae (6.11), 
the right-propagating wave g” changes into an exponential decreasing to the right and the left-propagating 
gf n wave changes into an exponential decreasing to the left, a fact that simplifies the identification of the 
right and left traveling and decaying waves in the actual implementation of the boundary conditions at z = 0 
and 2 = z max . 

It still remains to consider the case a > 16/3. For the positive values of fc 2 , we can introduce the 
wavelength A c = 27 r/fc c and for this range of a obtain A Jh z < y/3n/2. Thus, we see that a > 16/3 implies 


15 



a poor “points per wavelength” resolution even for the long waves A c > Ao = 27r/fco- This makes the 
choice a > 16/3 inappropriate for practical computations. Finally, regarding the last case that has not been 
considered yet, a = 0, we note that for this value of a equation (6.10) will have a double root q\ = </_ \ = 1. 
However, formulae (5.9) and (5.11) show that the case a = 0 4=>- k \ = 0 can be easily avoided by slightly 
changing the parameters of the discretization. 

For j — 2, we find from equation (6.10) that 




(6.13) 


Clearly, \q 2 \ < 1, \q 2 1 | > 1 for all relevant values of a (a > —3), i.e., the two components q 2 and q 2 n of 
(6.13) always correspond to evanescent waves. ■ 


6.2. Discrete one-way Helmholtz equations. In analogy with the continuous description in Sec- 
tion 4.2.1, we now construct the discrete one-way Helmholtz equations based on the solution (6.3) of the 
homogeneous finite-difference scheme (6.2). From the very beginning, we think of these discrete one-way 
Helmholtz equations as the relations to be used as boundary conditions for equation (6.1). 

According to Proposition 6.1, the discrete homogeneous equation (6.2) has four linearly independent 
eigenfunctions, two of which are either traveling or evanescent waves and two others are always evanescent 
waves; the presence of the latter (in contrast to the continuous case) is due to the fact that (6.2) is a 
fourth order finite-difference equation that approximates the original second-order differential equation. 
When constructing the discrete one-way Helmholtz equations, we, of course, first need to make sure that 
they handle the first pair of discrete waves, g" and in the same way that equations (4.9) handle 
the corresponding continuous waves. In addition, we need to decide how the discrete one-way Helmholtz 
equations will handle the second pair of discrete waves, q 2 and q 2 n , which are purely numerical (i.e., due 
to the use of a forth-order difference scheme). It is natural to require that the one-way-to-the-right discrete 
Helmholtz equation admit the right traveling/evanescent wave g” and the right evanescent wave q 2 and that 
the other two waves from representation (6.3), </f n (left traveling/evanescent) and q 2 n (left evanescent) be 
suppressed by this equation. Indeed, q^ n may either be traveling “the wrong way” or grow without bound 
as n — > +00 and q 2 n will always grow without bound as n — > +oo . 6 Clearly, if we use the one-way-to-the- 
right equation that possesses such properties as boundary condition for (6.2) near n = N, it will guarantee 
that the corresponding far-field solution ( n > N ) always be bounded at infinity and also that this solution 
may only be composed of outgoing (right propagating and/or right decaying) waves. In other words, the 
one-way-to-the-right discrete Helmholtz equation implies that in the far field n > N one can represent the 
solution u n in the “restricted” form 


u n = ciq? + c 2 q 2 , 


(6.14) 


as opposed to the general form (6.3). Formula (6.14) is equivalent to requiring that the vector 
[un- 3 ,un- 2 ,Un-i,un] be a linear combination of the two vectors [l,qi,ql,qf] and [1 ,q 2 ,q 2 ,q 2 ], which 
is the same as requiring that 


UjV-3 un- 2 Un- i u N 

Rank 1 q 1 qf qf = 2 . (6.15) 

. 1 Q 2 ql q 2 . 


6 Besides being “natural,” this choice is also motivated by the well-posedness considerations, as the analysis of [13, 20] 
suggests. 
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Relation (6.15) immediately yields the following two linearly independent conditions 



UN-3 

UN-2 

un- 1 


UN-2 

un- 1 

un 

det 

1 

qi 

ql 

= 0 , det 

'll 

qf 

ql 


1 

<?2 

ql . 


. <?2 

qi 

ql . 


which reduce to 

qiq 2 u N - 3 - (qx + q 2 )u N - 2 + u N - 1 = 0 (6.16a) 

and 

qiq 2 u N - 2 - (qi + q 2 )u N - 1 +u N = 0 . (6.16b) 

The two scalar equations (6.16a) and (6.16b) constitute the one-way-to-the-right discrete Helmholtz equation. 

The one-way-to-the-left discrete Helmholtz equation is constructed similarly. Symmetrically to the 
previous case, we require that it admit the left traveling/evanescent wave q and the left evanescent wave 
q^ n and that the other two waves from representation (6.3), q " (right traveling/evanescent) and q% (right 
evanescent) be prohibited by this equation. (From the standpoint of physics the two waves, q± n and q^ n , 
account for the phenomenon of backscattering.) The waves q ” and q% are to be suppressed in this case 
because q™ may either be traveling “the wrong way,” i.e., to the right, or grow without bound as n — > — oo 
and #2 will always grow without bound as n — l -oo. If the one-way-to-the-left discrete Helmholtz equation 
is used as boundary condition for (6.2) near n = 0, it will guarantee that the corresponding far-held solution 
(n < 0) always be bounded as z -¥ -oo, and also that this solution may only be composed of outgoing 
(left propagating and/or left decaying) waves. In other words, the one-way-to-the-left discrete Helmholtz 
equation implies that in the far held n < 0 one can represent the solution u n in the “restricted” form 

u n = c-iq[ n n + c- 2 q^ n , (6.17) 

as opposed to the general form (6.3). To make sure that representation (6.17) hold, we require that the 
vector [wo,Ui,U2, W3] be a linear combination of [l,^ 1 ,^ 2 ,^ 3 ] and [l,^ 1 ,^ 2 ,^ 3 ]: 

16 0 Ui u 2 u 3 

Rank 1 q ~ 1 q~ 2 q~ 3 = 2 . (6.18) 

. 1 q 2 l q 2 2 q 2 3 . 

Relation (6.18) is equivalent to the following two linearly independent homogeneous conditions: 

uo ~ (qi + q 2 )ui + qiq 2 u 2 = 0 (6.19a) 

and 

Ui - (qi + Q 2 )u 2 + qiq 2 u 3 = 0 , (6.19b) 

which constitute the one-way-to-the-left discrete Helmholtz equation. 

We note that splitting the general solution (6.3) into right- and left-going waves (equations (6.14) and 
(6.17), respectively), and allowing for only one direction while prohibiting the other at the corresponding 
edges of the interval constitutes the radiation principle in the finite-difference discrete framework. 
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Having constructed the one-way discrete Helmholtz equations (6.16) and (6.19), we now implement 
them as boundary conditions for the discrete homogeneous equation (6.2). If we consider the finite grid 
n = 0, 1, . . . ,N on the interval [0,2: max ], the five-node difference stencil cannot be centered at the near-edge 
nodes n = 0, 1, N — 1, and N. As a consequence, the number of equations in the linear system is less than 
the number of unknowns by four. To make the number of equations and the number of unknowns equal, 
we supplement equations (6.2) on the grid n = 2,3, . . . , N — 2 by equations (6.19a) and (6.19b) for n = 0 
and n = 1, respectively, and by equations (6.16a) and (6.16b) for n = A r — 1 and n = N, respectively. In 
doing so, we arrive at the following linear homogeneous algebraic system with N + 1 equations and N + 1 
unknowns: 






Au = 

0 , 




(6.20) 

where 


1 

~(Qi +Q2) 

6162 

0 

0 


0 




0 

1 

-(91 +62) 

9l92 

0 


0 




-1 

16 

(12 a - 30) 

16 

-1 


0 


A = 

1 

12 hi 

0 


-1 

16 

(12a - 30) 

16 

-1 

(6.21) 



0 


0 

9i9 2 

-(?i + 62) 

1 

0 




0 


0 

0 

9l92 

-(9i + 92 ) 

1 


and, obviously, u 


,u n ] t . 








The following Proposition 6.2 establishes the solvability and uniqueness of the solution for the non- 
homogeneous counterpart of system (6.20). 


Proposition 6.2. The linear non-homogeneous system of equations Au = f with the matrix A given 
by (6.21) is uniquely solvable for any right-hand side f = [/o, fi , • • • , /jv] t • 

Proof. We show that the corresponding homogeneous system (6.20) has only trivial solution. Indeed, the 
only solution to any of the equations of Au = 0 except the first two and the last two is a linear combination 
of the type (6.3). However, each of the components of (6.3) is explicitly eliminated by one of the boundary 
conditions (6.16a), (6.16b), (6.19a), or (6.19b), i.e., by one of the one-way discrete Helmholtz equations (the 
first two and the last two equations of Au = 0). Therefore, the only solution to the homogeneous system is 
the trivial one. 7 ■ 

Although we have just shown that one can find the solution to Au = /, for any given f = [/o, /i, • • • , /tv], 
this solution will not, in fact, correctly approximate the corresponding solution of the non-homogeneous 
differential equation, or in other words, will not, generally speaking, be the discrete radiation solution from 
the sources / = [/o,/i,--- ,/jv]- The reason for this discrepancy is that the one-way Helmholtz equations 
which are used in the first two and last two rows of the matrix A have been constructed for the homogeneous 
case. As a result, these four equations will not handle correctly the near-boundary source terms, which may, 
generally speaking, be present. The “cure” to this problem, in the form of a a local modification to /, is 
derived in Section 6.4. 

7 This solvability result is obviously similar to the one in the continuous case, see Section 4.2.1. 
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In our simulations (see Section 7), we solve, the finite-difference Helmholtz equation by inverting the 
matrix A of (6.21). However, for the purpose of deriving the two-way ABCs that would correctly handle the 
near-boundary inhomogeneities, we now show how to construct the solution u by using the Green’s function 
of the finite-difference operator of (6.2). As we shall see, this approach is rather expensive numerically and 
thus not useful for actual computing. However, it provides the most conceptually straightforward way to build 
the radiation solution. Moreover, the analysis that employs the Green’s function reveals the mechanism of 
the aforementioned discrepancy between the radiation from the sources f = [/o, /i, ■ ■ ■ , /.v] and the solution 
to Au = /. 


6.3. Radiation solution by means of the Green’s function. In this section, we introduce a 
problem very similar to (6.1), except that the solution u is now defined on the infinite grid n = 0, ±1, ±2, . . . , 
and the right-hand side f n is compactly supported: 


Un — 2 T 16u n _i T 16'U n _|_i W n _|_ 2 . 2 


12 hi 

fn — 0 


for 


+ — fn •> ^ — 0, il, i2, ... , 


n < 0 and n > N . 


(6.22) 


We also require that the solution u n of (6.22) satisfy the radiation principle in the areas of homogeneity 
n < 0 and n > N. In other words, we require that for n < 0 one can represent u n in the form (6.17) and for 
n > N in the form (6.14). This is the most general formulation of the problem of finding the solution that 
corresponds to the radiation of waves by the sources / = [/ 0 , /i, . . . , /jv] T in the finite-difference framework. 

To solve this problem, we introduce the fundamental solution G n (free-space Green’s function) for the 
one-dimensional discrete Helmholtz operator, which is defined on the entire infinite grid n = 0,±1,±2, ... 
and is the solution of the equation 

_ G n-i + i6G n- i _ 3QG n + l6G n+1 - G n+2 


12 hi 


+ k 2 c G n = 5 n , n = 0, ±1, ±2, . . . , 


(6.23) 


where 


S n = 


1, n = 0 
0, n ^ 0 


We also require that the Green’s function G n satisfy the radiation principle as n — > ±oo, or in other words, 
that it can be represented in the following form: 


G n = 


«i l l\ + n>0 


l Mr n + n < ° 

Proposition 6.3. The values of the constants a±, a?, b\, b 2 in (6.24) are given by 

12h 2 q\ 


ai = 


(? 2 1 - tfiXtfi 1 - 9 i )(«2 - qi ) 


(6.24) 


(6.25a) 


02 


— 12h 2 z q2 

{qf 1 - « 2 )(«r 1 - «2)(<72 - qi) 


(6.25b) 


bi 


-V2h 2 q-' 

feT 1 - qf v )(qf' - ftOter 1 - qi ) ’ 


(6.25c) 
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(6.25d) 


i _ 1^ h 2 z q 2 

(Q 2 1 - <h l )(<h' - Q2)(q2 1 - 4l) ’ 

Proof. To find these four constants, we need four equations. By matching the two branches (6.24) of the 
Green’s function G n at n = 0 we immediately obtain one equation 

ui + 12 = bi + 62 • (6.26a) 

The other three equations for the coefficients of (6.24) are obtained from the original equation (6.23) written 
for the nodes n = 0, 1 and —1. For n = 0we have 

-G~ 2 + 16G -1 + (12c* - 30)G'° + 16G' 1 - G 2 — 12 h\ , 


or 

— (biQi + &24 2 ) + 16 (biqi + 6242) + (12c* — 30)(oi + 02) + 16(fliQi + 02Q2) — ( a iq 2 + 0242) = 12h 2 • 

The previous equation can be simplified by subtracting from it the following relation 

— (oi4i 2 + °2<?2 2 ) + 16 (oi4i 1 + I 2 Q 2 *) T (12c* — 30) (01 + 02 ) + 16(oiQi + 0242 ) — (oi4i + = 0 , 

which comes from the fact that each branch of the Green’s function (the right branch dig” + 02^2 i n this 
particular instance) satisfies the homogeneous finite-difference equation (6.2). The subtraction yields: 

— (Mi + M2) + 16 (Mi + M2) — 16 (ai<?x 1 + 02^2 1 ) + (°i4i 2 + a 2 Q 2 2 ) = 12h 2 • (6.26b) 

For n = 1 equation (6.23) takes the form 

-G -1 + 16G° + (12a - 30)G J + 16G 2 - G 3 = 0 

and again, using the homogeneous equation for the right branch of the Green’s function, we obtain 

- (biqi + M2) + {aiqf 1 + a^q^ 1 ) = 0. (6.26c) 

Finally, for n = —1 we have 

-G“ 3 + 16G -2 + (12a - 30)G _1 + 16G° - G 1 = 0 . 

Combining this relation with the homogeneous difference equation for the left branch of the Green’s function, 
we arrive at 


(biqi 1 + M 2 *) - (ai9i + «'2<?2) = 0 . (6.26d) 

Now we need to solve equations (6.26) for oq, 02, 61, 62- First, we multiply (6.26c) by 16 and substitute 
it into (6.26b) and then rewrite all four equations as follows 


<?1 2 42 2 -Ql -qI 


O-i 


1 

I— 1 

to 

G lO 

1 

Qi 1 Q 2 1 -Qi 0.2 


0.2 


0 

11-1-1 


bi 


0 

4i 42 -4i _1 -42" 1 . 


. b 2 . 


0 
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The determinant of system (6.27) is easily reduced to a Vandermonde determinant, which eventually leads 
to expressions (6.25). ■ 

From the definition of G”, we have the following 

Proposition 6.4. For any given right-hand side /„ compactly supported on [0,1 ,7V], the solution 

to (6.22), subject to the radiation principle, is given by the convolution 

m=N 

u n =J2 fmG n ~ m , n = 0, ±1, ±2 (6.28) 

m = 0 


6.4. Radiation solution by means of inverting the matrix A. The cost of calculating the con- 
volutions in (6.28) for n = 0, 1, . . . , N is 0(N 2 ). We now show that the portion of the solution (6.28) that 
we are interested in, namely, u n for n = 0, 1, . . . , N , can be recovered by means of inverting the matrix A 
of (6.21). The cost of this inversion will be only O(N) operations because the matrix A is pentadiagonal, 
see Section 6.7 for additional detail. 

Proposition 6.5. Let A be defined by (6.21) and u = [«o,«i,... ,mjv] t be defined by (6.28) for 
n = 0, 1 , . . . ,N. Denote f = [A, A, • • • > Av-i, Av] T - Then, Au = f , where 


0 


fo 

0 


A 

A 


0 


+ 


f N—2 


0 

0 


Av-i 

0 


In 


(6.29) 


f def 1 

/0 “ mf 


(f 0 G° + fiG' 1 ) - ( qi + q 2 ) (/oG 1 + fiG°) + qi q 2 (/ 0 G 2 + AG 1 ) , 


(6.30a) 


f def 1 

h ~ 12h f 


(/oG 1 + AG° + AG" 1 ) - (q, + q 2 ) (AG 2 + A G 1 + f 2 G°) + 


qiqi (AG 3 + AG 2 + AG 1 ) , 


(6.30b) 


r def 1 

JN-1 — 


12 h 2 . 


QiQz ( Av— 2G 1 + Av-i G 2 + A vG 3 ) — 


(^i + qi ) (AV- 2 G 0 + /jv-iG 1 + AvG 2 ) + (AV- 2 G 1 + /jv-i G° + AvG 1 ) 


(6.30c) 


and 


r def 1 

Av = 


12 hi 


q\ q 2 (/n-iG 1 + AvG 2 ) — ( qi + q 2 ) (Av-iG° + AvG 1 ) + (Av-iG 1 + AvG 0 ) 


(6.30d) 
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Proof. By definition of the Green’s function G n (see Section 6.3), (Au) n — f n for 2 < n < N — 2. 
Indeed, for 2 < n < N — 2 we have 

12 h z — Vjji—2 “h 16tx n _i “I - (12 ol 30)u n ^n -\- 2 

N N N 

= - E fmG n - 2 - m + 16 E fmG n - x - m + (12a - 30) E fmG n ~ m 

m = 0 m = 0 m = 0 

N N 

+ 16 E fmG n+1 ~ m - E fmG n+2 ~ m 

m = 0 m = 0 

N 

= E /"» (— G n_2_m + leG”- 1 -™ + (12a - 30)G ra_m + l&G n+1 ~ m - G n+2 - m ) 

m = 0 

JV JV 

= E 12 h 2 z 5 n - m fm = 12 h 2 z E fm^n—m = l2h 2 J n . 

m = 0 m=0 

As for (4 m) 0 , ( j 4m) 1 , (Au) n _ 1 , and ( Au ) N , these four components need to be calculated separately. 
They will, generally speaking, differ from /o, /i, /jv-i, and /jv, respectively, because of the special structure 
of the first two and the last two rows of the matrix A , which admit waves going in only one direction, see 
Section 6.2. 

We start the analysis from the left edge of the interval Clearly, any f m for m > 2 is not going to 
contribute to (1«) 0 because when substituting u of (6.28) into (6.19a) we, in fact, substitute only the left 
branch of the Green’s function G n ~ m , see (6.24). Indeed, in formula (6.19a) we only need the values of u n 
for n = 0, 1, 2, and if m > 2 this implies n — m < 0. The left branch of the Green’s function (6.24) by 
definition turns (6.19a) into an identity, therefore ( Au) 0 is not affected by f m for m > 2. Consequently, 

(Au) 0 = {A [f 0 G n + fiG n ~ 1 ]) 0 , 

which proves (6.30a). Similarly, substitution of the left branch of the Green’s function into (6.19b) suggests 
that any f m for to > 3 is not going to contribute to {Au) 1 . Therefore, 

{Au\ = {A [foG n + /iG” -1 + hG n ~ 2 ]) 1 , 

which proves (6.30b). 

Similar analysis is conducted for the right edge of the interval. Only /jy and /jv-i affect (Au)n = /v 
because for all other components of the R.HS / the contribution to the solution u at n — N — 2, N — 1, N is 
given by the right branch of the Green’s function only; then the explicit form of the solution (6.28) and the 
definition of A (6.21) easily yield expression (6.30d). Analogously, only three components of the right-hand 
side, /jv, f N—i , and /jv-2, contribute to (Am)jv-i = /n- i, which together with (6.28) and (6.21) implies 
(6.30c). ■ 

From the standpoint of the original physical model the situation near z = z max differs substantially 
from the situation near z = 0, because we can always make the effect of nonlinearity and/or variation of 
coefficients near 2 = 2 max negligible, by taking 2 max sufficiently large. Therefore, from here on we will always 
assume that /at = fiy-i = fN - 2 = 0. Obviously, if we use the RHS / = [/o , / 1 , • • • ,/jv-3,0,0,0] t of this 
particular kind as source terms in (6.22), then for the corresponding solution u = [uo,Ui, . . . ,ujv] we will 
have (Am) jv _ 1 = fN-i = 0, see (6.30c), and (Au) N = /jv = 0, see (6.30d). In other words, the modified 
hight-hand side / of (6.29) in this case becomes / = [/o, / 1 , / 2 , • • • ,/at_ 3,0,0,0] t . 
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Let us emphasize that /o = (Am) 0 , see (6.30a), depends on / 0 and /i, and /i = ( Au)i , see (6.30b), 
depends on / 0 , fi, and fi . Likewise, in order to obtain /jv-i = (^w)at_i = 0, see (6.30c) and /at = 
( Au) n = 0, see (6.30d), in addition to the obvious requirement that /at = Jn-i = 0, we also need to impose 
f N—2 — 0. 

Propositions 6.2 and 6.5 guarantee that the only solution of the linear system Au = /, where / = 
[/o, A, h, ■ ■ • , In- 3,0, 0, 0] T , is the solution u of (6.22) with the RHS / = [/ 0 , fuh,--- , In- 3,0, 0, 0] subject 
to the radiation principle. Thus, we have addressed the concern raised in the end of Section 6.2, namely, 
which modifications to the right-hand side / are needed so that the solution obtained by inverting the matrix 
A will coincide with the pure radiation solution from these particular sources /. Provided that near the 
right edge of the interval the RHS is zero: /a? = /n - i = /at - 2 = 0, it turns out that these modifications are 
local and require only the replacement of the two old quantities /o and /i near the left edge of the interval 
by the new quantities /o and /i, respectively. It is also important to mention that formulae (6.30a), (6.30b) 
are by themselves local as well, and therefore the overall modification / i — ? / amounts to only local, and 
thus numerically inexpensive, operations on the grid near n = 0. 

6 . 5 . Adding the incoming power. The boundary conditions at z = 0 should guarantee the complete 

transparency of this boundary for all backscattered waves and at the same time be capable of accurately 
prescribing the incoming signal; the combination of these two properties has been referred to as the two- 
way ABCs. Similarly to the continuous case analyzed in Section 4.2.2, the incoming signal u? nc results in 
a forward propagating wave, given by The grid function v n = u° nc qi solves all equations of the 

homogeneous system Av = 0 except for the first two, which are the one-way-to-the-left discrete Helmholtz 
equation (6.19). Therefore, by applying the matrix A of (6.21) to the vector v we create a right-hand side 
that we denote by g. It is easy to see that 

1 - (<?i +q 2 )qi +qfqz 
4i(l - (Qi +q 2 )qi +qfq 2 ) 

0 

0 

Proposition 6.2 guarantees that the only solution of the system of equations Av = g, where g is given by 
formula (6.31), is v = u® nc q™. Note, the inhomogeneity g of (6.31) is a discrete counterpart of the right- 
hand side of relation (4.12) (and (4.11)) obtained when introducing the incoming signal in the continuous 
framework, see Sections 4.2.2 and 4.2.3. 

6.6. Obtaining the overall solution. We can, finally, put together the foregoing stages of the deriva- 
tion. Assume that there is a given RHS of the original equation (6.1) / = [/cn fi, /2, • • • ,/Ar-3,0,0,0] T . 
To obtain the solution with the incoming power Vr^ K q\ l added, we first construct the new RHS / on the 
basis of / according to formulae (6.29) and (6.30a), (6.30b). Then, we construct the additional source terms 
g according to formula (6.31). Due to the linear superposition principle and according to Proposition 6.2 
that guarantees solvability and uniqueness, we immediately see that the grid function u that we recover by 
solving the overall system 

Au = f + g , (6.32) 

is the solution that we are looking for. Indeed, including / on the right-hand side of (6.32) guarantees the 
radiation from the original sources / both to the left and to the right and including g on the right-hand side 


(6.31) 
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of (6.32) guarantees that the correct incoming signal will be added. The system (6.32) is, of course, 

solved by inverting the matrix A only once and not by solving separately with the RHSs / and g. 

Thus, setting the desired boundary conditions at z = 0 and 0 = 3 ma x is reduced to building and inverting 
the special matrix A of (6.21) and also modifying the right-hand side of the equation / \ — > / + g. We again 
emphasize that the latter modification is not computationally expensive as both f and g are obtained by only 
local operations on the grid near n = 0. These operations will come at virtually no cost when implementing 
the algorithm numerically. 

To conclude this section we note that the solvability and well-posedness analysis of general one- 
dimensional systems of finite-difference equations can be found in [13,20]. 

6.7. Solution of Au = / + g. We solve the system Au = / + g using standard LU decomposition; 
for a pentadiagonal matrix A the components of this decomposition will obviously be banded as well. As 
the equation Au = f + g needs to be solved many times with changing source term but with the same A, 
at the beginning of a simulation we calculate once the LU decomposition of A, and use it throughout the 
iterations. Therefore, the costs per iteration in terms of solving this equation are only due to the backward 
substitution, which is O(N) arithmetic operations. 

7. Numerical experiments. To assess the numerical performance of our algorithm, we first solve a 
linear problem with variable coefficients in several different settings. 

7.1. Linear problem with variable coefficients and backscattering. On a slender rectangular 
domain in the (r, z) coordinates, [0, r max ] x [0, 2 max ], where r max = tt/2 is fixed, and 2 max will vary as an 
essential part of testing the methodology, we are recovering the following solution: 


A — -Uright A C • A] eft , 


(7.1) 


where C is a constant, and the right and left propagating components .Bright and Bi e f t are given by: 

B right = eV*o 2 -" 2 * cos(i'r) [l + ez 4 e~ z ] , (7.2a) 

Bi e ft = e -V*o-" 2 * cos (isr)e-( z/ ^ 2 . (7.2b) 

In the framework of our study, the left propagating component Bi e f t of (7.2b) is interpreted as backscattering. 
Several parameters that control the actual shape of the solution (7.1) are: k o is the wavenumber that 
corresponds to the homogeneous medium, see Sections 1 and 2; v is the transversal frequency; e in (7.2a) 
determines the extent of deviation from the constant-coefficient case for the right propagating mode (see 
below); and (3 in (7.2b) determines the spatial (longitudinal) extent, to which the backscattered waves are 
present in the solution. In the linear case we, of course, introduce the backscattered waves artificially, but 
we are trying to follow the physically interesting situation when these waves are generated inside the domain 
and propagate toward and through the left boundary z = 0. The constant C is introduced in (7.1) so as to 
control the magnitude of the backscattered signal relative to the forward propagating signal and in particular 
to be able to fully eliminate backscattering ( C = 0) if desired. 
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Substituting B r i g ht of (7.2a) into equation (4.5a), we obtain: 

Alright + ^o^right = ee l ^ k °~ l '~ z cos (vr)e~ z z 2 2 k q — v 2 (4 z — z 2 ) +12 — 8 z + z 2 


e Z z 2 


2i^/k^ — v 2 (4 2 — z 2 ) + 12 — 82: + . 


1 + e • z 4 e 


4 P — z 


L E ri 


(7.3) 


right 


— -Fright * Fright • 

We therefore conclude that B r i g ht of (7.2a) satisfies the variable-coefficient equation 

A -Bright 4“ ^'right( 2 ') -bright — 0 i 

where A' 2 ight (;2) = fcg(l + eB r i g ht(^)) and B r i g ht(;2) is defined by equalities (7.3). We indeed see that e controls 
the extent of spatial variation of the wavenumber A+ght- The solution B r i g ht is driven by the incoming wave 

B inc = e i v fk °~ 1 ' 2 z cos (yr) , * < 0 . (7.4) 


Similarly, the backscattered solution Eeft of (7.2b) satisfies the variable-coefficient equation 

A£) e f t + k 2 e{t (z)E\ e f t = 0 , 

where kf eft (z) = &§(! + Bi e ft(z)) and 


Ai 


Tleft(-) — i 2 

L 

For the overall solution E of (7.1) we obviously have 




;s_ _ 4 ^ 

> 2 /3 2 + 0 2 


AB + k 2 (z)E = 0 , 


(7.5) 

(7.6) 


where 


k 2 (z) = kl 


Bright 


right. 


2 C ■ Bi e ft 

Ki, - 


''left. " 


B 


The driving incoming signal for equation (7.6) is B inc of (7.4), evaluated at z = 0. The variable-coefficient 
linear equation (7.6) for B will be solved on the domain [0,r max ] x [0,2 max ] with the homogeneous radiation 
boundary condition (4.10a) at z = z max and non-homogeneous (two-way) radiation boundary condition (4.12) 
at 2 = 0. The boundary conditions at r = 0 and r = r max are symmetry and zero Dirichlet, respectively, 
which corresponds to the general construction of Section 5, as well as the particular explicit form of the 
solution (7.1), (7.2) that we use here. The solution will be obtained by iterations described in Section 4.2; 
the corresponding discrete solution methodology is delineated in Sections 5 and 6. 

Our primary goal when solving numerically the foregoing linear problem is to demonstrate that the 
algorithm that we have constructed indeed possesses the design properties, i.e., (1) converges with the 
fourth order of accuracy when the grid is refined, and (2) properly handles the radiation of waves (including 
backscattering) or in other words, introduces no reflection from the boundaries z = 0 and 2 = z max back into 
the domain. A secondary goal is deriving the guidelines for subsequent nonlinear simulations, for example, 
how geometric parameters, such as domain size, may affect the solution. 

The forthcoming series of computational experiments corroborates our expectations in terms of grid 
convergence and handling the backscattered waves, and also provides for a comparison between the following 
two algorithms: The one constructed in this paper with the two-way ABC at the boundary z = 0 and a 
more traditional one with the Dirichlet boundary condition at z = 0 (at the far-field boundary z = 2 max we 
set the same radiation ABC in both cases). 
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7 . 1 . 1 . Traditional approach — Dirichlet boundary condition. The algorithm that we have just 
referred to as a more traditional one is formulated with the Dirichlet boundary condition for E at z = 0. In 
fact, already from the standpoint of physics one can anticipate that this algorithm is not going to perform 
well when backscattering is present. Indeed, the physical setup of the model implies that all the information 
available at 2 = 0 pertains only to the incoming wave. Thus, we basically cannot say anything about the 
backscattered signal ahead of time because it is generated inside the domain (in the current example we, of 
course, know everything because we simply construct a sample solution including the backscattering, then 
produce the corresponding sources/inhomogeneities, and finally recover the same solution by the numerical 
method, but this is done only for the demonstration purposes.) When constructing the two-way ABCs, 
we do not make and do not need any assumptions regarding the backscattered wave, we simply make the 
boundary transparent for all such waves. In contrast, in the Dirichlet case we can only specify the incoming 
wave as the boundary data because no explicit information about other waves is available. Mathematically, 
this amounts to making the following assumption/approximation: 

E(r,0) = E? nc (r) , (7.7) 

which, in contradistinction to (2.1), prescribes the entire field at z = 0, rather than its incoming component 
only. Consequently, the Dirichlet boundary condition will essentially reflect all backscattered waves reaching 
2 = 0 back into the medium, in contrast with the two-way ABC, which will let them go through. We 
therefore expect that the algorithm with the Dirichlet boundary condition (7.7) at z = 0 may produce 
reasonable results only if no backscattered waves are present in the solution. Otherwise, the error should 
be roughly of the magnitude of the backscattered signal. The numerical results below corroborate these 
expectations. 

Note that to enforce the Dirichlet boundary condition at z = 0 for the discretization we obviously assign 
a prescribed value to the solution at the leftmost grid node n = 0. Besides, in the framework of the fourth- 
order scheme that we are using, we need an additional relation to be specified right next to the boundary 
at n = 1. This is similar to obtaining the discrete one-way Helmholtz equations in the form of two scalar 
relations , see Section 6.2. The additional relation for the Dirichlet boundary conditions should be merely an 
approximation of the underlying differential equation at n = 1, but this cannot be the same approximation 
that we are using for the interior nodes (n > 2) because the latter employs a five-node wide symmetric stencil. 
Thus, either a one-sided difference or a compact Pade-type approximation needs to be used at n = 1. We 
chose the fourth-order Pade [6] on a three-node wide stencil in the particular form proposed in [25] because 
as opposed to the “long” non-symmetric differences, it preserves the pentadiagonal structure of the matrix. 

7 . 1 . 2 . Results. For the simulations in the linear case we have chosen the following particular values 

of parameters (see formulae (7.1), (7.2)): = 20, e = 0.2, v = 3 or v = 1, z max = 30 or z max = 10, 

(3 = 3, C = 1/2 for the case with backscattering, and C = 0 for for the case with no backscattering. The 
wavelengths in the r and z directions are A r = 2-it/v and \ z = 2n/ko , respectively. We choose the grid sizes 
h r and h z accordingly as fractions of the corresponding wavelengths: For the grid convergence study we 
refine the grid synchronously in both r and z directions. We note that having the same resolution (nodes per 
wavelength) in both directions yields the cell aspect ratio of h r /h z = A r /A z = ko/v , which in our simulations 
is equal to either 20/1 or 20/3. 

We have looked at the values of the relative error (the difference between the computed and exact 
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solution normalized by the maximum of the exact solution over the domain) in the maximum norm: 


Error = 


max | ©computed ©exact | 

(r>z) 

max I -©exact | 

(r,z) 


(7.8) 


The results are summarized in Tables 7.1 and 7.2 for v = 1 and v = 3, respectively. In both tables all values, 
except those in the last column, correspond to 2 max = 30. 


Table 7.1 

Maximum relative error (7.8) of the calculated solution in the linear case for v — 1 . 


Grid sizes 

Backscattering 

Off (C = 0) 

On (G = l/2) 

Boundary condition at z = 0 

Dirichlet 

Two-way 

Dirichlet 

Two-way 

Two-way 

^max — 30 

^max — 10 

h r = A r /10, h z = A z /10 

0.256 

0.257 

0.33 

0.24 

0.16 

h r = A r /20, h z = A*/20 

0.0165 

0.0165 

0.33 

0.016 

0.01 

h r = A r /40, h z = Aj/40 

0.001 

0.001 

0.33 

0.001 

0.0012 

h r = A r /80, h z = A*/80 

6.5- 10- 5 

6.5- 10- 5 

0.33 

6.5- 10~ 5 

0.00075 


Table 7.2 

Same as Table 7.1, with v = 3. 


Grid sizes 

Backscattering 

Off (C = 0) 

On (£7=1/2) 

Boundary condition at z = 0 

Dirichlet 

Two-way 

Dirichlet 

Two-way 

Two-way 

^max — 30 

^max = 10 

h r = A r /10, h z = Az/10 

0.25 

0.25 

0.33 

0.24 

0.089 

h T = A r /20, h z = A s /20 

0.016 

0.016 

0.33 

0.015 

0.0064 

h T = A r /40, h z = Aj/40 

0.001 

0.001 

0.33 

0.001 

0.0012 

h r = A r /80, h z = Aj/80 

6.3 ■ 10~ 5 

6.3 • 10“ 5 

0.33 

6.3- 10~ 5 

0.00075 


From Tables 7.1 and 7.2 we first conclude that, as expected, the Dirichlet boundary condition (7.7) 
provides no convergence when the backscattering is present (third column). In all other columns we observe 
a fourth-order grid convergence, because every time the grid is refined by a factor of two in each direction, 
the value of the error drops by approximately a factor of sixteen (except for the last column of each table, 
which will be discussed later). Thus, the algorithm that we have constructed indeed possesses the design 
convergence properties. Besides, we clearly see that the left propagating waves in the solution present no 
problem from the standpoint of numerics for the algorithm with the two-way ABC at z = 0. 

Let us now return to the data appearing in the rightmost columns of both Table 7.1 and Table 7.2. 
These data clearly do not demonstrate the fourth-order grid convergence. The only difference between these 
data and all other data in the tables is that the rightmost columns correspond to a smaller computational 
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domain in the z direction, £ max = 10, as opposed to z max = 30. The breakdown of the grid convergence that 
we observe on the small domain has the following explanation. 

The boundary condition that we specify at z = z m ax is the homogeneous radiation boundary condition 
(4.10a), which is approximated by the one-way-to-the-right discrete Helmholtz equation (6.16). Both the 
continuous (4.10a) and discrete (6.16) radiation boundary conditions at z = z max were obtained under the 
key assumption that the governing equation near z = 2 max reduce to the constant-coefficient Helmholtz 
equation A E + kgE = 0. In other words, this means that the mode flight of (7.2a) has to reduce to 
the “pure” propagating mode e*''/*? - " 2 * cos (z/r), and the mode -E) e ft of (7.2b) has to effectively vanish at 
z — 2 m ax- From the specific form of the modes that we analyze, see (7.2), we conclude that the larger we 
take the domain [0,z max ] the better the quality of the agreement with the desired properties near 2 = 3 max 
is going to be. In other words, for the smaller domain 2 max = 10 we are essentially trying to apply a 
homogeneous radiation boundary condition to the equation, which is not “sufficiently homogeneous” itself 
and therefore, the error is dominated by this discrepancy, rather than the actual truncation error associated 
with the discretization of the differential operator. As a consequence, we do not observe the fourth-order 
grid convergence for the smaller domain. This demonstrates the importance of choosing z max sufficiently 
large, so that the homogeneous radiation boundary conditions can be applied successfully. 



Fig. 7.1. Behavior of the error (7.9) for v = 1, two-way ABC at z = 0, h r = A r /20, h z = Az/20, f3 = 3, and z m ax = 30. 


Another interesting phenomenon that we would like to discuss in the framework of the linear case is the 
behavior of the error as a function of the coordinate z. A typical example in Figure 7.1(a), which corresponds 
to the case of no backscattering, shows a linear growth of the error with z except in the area of a small 
“bump” near the boundary z = 0. The actual quantity represented in Figure 7.1(a) is 


Error(;?) 


max | E computed E ( 


exact I 


max | inexact | 


(7.9) 


A similar error pattern is obtained for the case with backscattering, as shown in Figure 7.1(b). The curve 
in Figure 7.1(b) can be described as an oscillatory region next to the boundary z = 0 associated with 
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backscattering (the magnitude the error is still small there) followed again by a stretch of linear growth. 

It is, in fact, easy to see where this linear growth comes from. Proposition 6.1 implies that the discrete 
right propagating mode approximates the continuous right propagating mode e zkcZ = e tkc ' hzn (in the 
notations of this section, k % = ^A’o — p2 )- Indeed, assuming that k c ■ h z is small, we have obtained that 
qi = e lkc ' hz + O ((fc c • h z ) 5 ), see formula (6.4a). Consequently, under the same assumption we have g” = 
e ik c -h z n _|_ q (n(k c • h z ) 5 ) = e lkcZ + O ( zh 4 ) because z = h z n. As 0 < z < 2 max , we see that the error grows 
linearly in 2 and that the maximal error is 0(z m & x ■ h 4 ). The aforementioned linear growth of the error 
explains, in particular, why on coarser grids we obtain smaller maximal error for 3 max = 10 (fifth column) 
than for z max = 30 (fourth column), see Tables 7.1 and 7.2. 

It is, in fact, instructive to see how the 
error curve similar to those displayed in Fig- 
ure 7.1 would look for a solution computed on 
the small domain z max = 10. In Figure 7.2 
we show such a curve for exactly the same 
set of parameters used for computations that 
led to Figure 7.1(a), except that 2 max is equal 
to 10 instead of 30. Although the magnitude 
of the error is small, we observe oscillations 
throughout the entire domain. As we have no 
backscattering in this case (C = 0), the os- 
cillations may come only from the right (far- 
field) boundary z = 2 max . In fact, these oscil- 
lations are an early manifestation of the phe- 
nomenon that we have discussed earlier. On 
small domains, the application of the homoge- 
neous far-field radiation boundary conditions 
(4.10a.) and (6.16) is not fully “legitimate” be- 
cause the governing equation itself is not suf- 
ficiently close yet to the constant coefficient version A E + k^E = 0. The inconsistency gives rise to the 
oscillations shown in Figure 7.2. For finer grids this inconsistency, as we have seen, prevents the methodol- 
ogy from converging on small domains with the theoretically prescribed rate of 0{h 4 ). 

7.2. Nonlinear problem. Having corroborated the design properties of the numerical algorithm in 
the linear regime in Section 7.1, we now address its performance for the nonlinear case. In all cases that we 
analyze hereafter we take the value of ko = 8 and as before denote A z = 2w /fc 0 . In addition, in all simulations 
the solution is driven by the incoming signal 

El c {r)= e ~ T \ (7.10) 

The key quantity in the NLS model, as far as nonlinear self-focusing and singularity formation are 
concerned, is the ratio of the power of E® nc and the critical power N c (see Section 3.1). Therefore, we now 
briefly review the calculation of the critical power for the NLS (3.3). 

7.2.1. Critical power. Weinstein [30] had proved that the critical power for singularity formation in 
the critical NLS, N c , is equal to the power of the so-called waveguide solution. In the case of the (l-l-l)D 


x 10 -3 



Fig. 7.2. Same as Figure 7.1(a) with z mAX = 10. 
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critical NLS (3.3), the waveguide solutions are of the form 

tp(z,r) = exp (iaz)Q(r; a) . 

Substitution of this solution in (3.3) shows that the waveguide profile Q satisfies: 

Q rr -aQ + Q 5 = 0 , Q'(0) = 0 , Q(oo) = 0 . 


Integration of this equation yields: 

Q(r;a) = (3a) 1 / 4 sech 1 / 2 (2 v / ar) . 

Therefore, a necessary condition for singularity formation in (3.3) is that 

/•OO 

/ \Mf)\ldf>N c , 

Jo 

where 


N r 


-f 


Q' 2 (r) dr = 




In dimensional variables, this condition is 


f 


\EL(r)\ldr> 


N c 

k 0 \/e ' 


Therefore, the fractional critical power of E® nc of (7.10) is 


. JfWdr [2 

Nc/ koy/t 


koV~e • 


(7.11) 



Fig. 7.3. Grid convergence for e = 0.04, z m ax = 20, r max /z max = 
1, h z = Aj/lO, for h r = A z /2 (solid line), h T = A z /4 (dotted line), 
and h r = A z /8 (dashed line). 


7.2.2. Results. We start with a moder- 
ate nonlinearity in equation (3.2), e = 0.04, 
which, according to (7.11), corresponds to 74% 
of the critical power when fc 0 = 8. Our goal 
is to first demonstrate the grid convergence of 
the algorithm. We also compare the two-way 
ABC against the standard Dirichlet boundary 
condition at 2 = 0, as we did in the linear 
case, both from the standpoint of accuracy of 
the solution and the rate of convergence of our 
iterative scheme. 

For the grid convergence study we first 
choose the following parameters: z max = 20, 
^max/^max — 1? hz ~ hr ~ -^z/2. In OUr 

computations we have observed that changing 
the discretization parameters in the r direc- 
tion may exert a more noticeable influence on 
the solution than changing the discretization 


in the z direction. Therefore, we initially refine the grid in the r direction only and in Figure 7.3 present 
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three solution curves: Absolute value |-B C omputed(0, z)\ on the axis of symmetry r = 0 as a function of 2 
for h r = \ z /2, X z /4, and X z /8. We see that the last two curves that correspond to h r = X z /4 and X z /8 
are virtually indistinguishable from one another and both differ noticeably from the first one obtained on 
a coarser grid h r = X z /2. We therefore conclude that as the grid is refined the numerical solution does 
converge, even so in this nonlinear case we do not know what the exact solution is and consequently cannot 
explicitly find the error. 

We note that we plot the values of the computed solution on the axis of symmetry r = 0 because this is 
the most interesting location in the domain where the genuinely nonlinear phenomena take place. A clear 
manifestation of this nonlinear phenomena is the “bump”, or peak, on the solution curve in Figure 7.3, whose 
value is higher than that of the incoming wave £-^(0) = 1. Clearly, in the absence of nonlinear effects (i.e., 
e = 0), an unfocused input beam, such as (7.10), would simply diffract while propagating to the right, i.e., 
toward large z's, with its maximum amplitude monotonically decreasing. The amplification of the incoming 
signal due to the nonlinear response of the medium is called self-focusing , and is well-known within the NLS 
framework. 




(a) hz — A 2 /IO, h r — X z /4 


(b) h z = A*/20, h r = Xz/8 


Fig. 7.4. Backscattering for e = 0.04, z max = 20 and r max /z max = 1. 


Another interesting phenomenon, which is actually the one that our methodology has been specifically 
designed to capture, is backscattering. In the previous linear studies in Section 7.1, the extent of backscatter- 
ing was predetermined by the value of C. To estimate the extent of backscattering in the current nonlinear 
case, we plot the quantity |£ C om P uted(A 0) — £?, c (r)| as a function of r. In Figure 7.4(a) we show the corre- 
sponding graph for e = 0.04, 2 max = 20, r max /z max = 1, h z = A z /10, and h r = X z /4. From Figure 7.4(a) 
we conclude that most backscattering occurs around the axis of symmetry r = 0, and that the magnitude 
of backscattering there is about 1.2% of the incoming power. Backscattering obviously accounts for the 
deviation of the solution curve at z = 0 in Figure 7.3 from the incoming signal value there, which is equal 
to 1. 

A comprehensive grid refinement study should, of course, include refinement in the 2 : direction along 
with the refinement in the r direction. In addition to the cases reported previously, we have run several 
others, refining the grid either separately in each direction or synchronously in both directions, and also 
changing the size of the computational domain. Note that determining the correct, i.e., sufficiently large, 
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size of the computational domain is important, because choosing it too small in the ^ direction may cause 
reflections from the boundary z = z max (Section 7.1), and choosing the domain too small in the r direction 
is dangerous because the boundary r = r max is reflecting and the reflections may, in fact, completely destroy 
the solution (we have actually observed the latter phenomenon in our computations). 

Basically, the solutions that we have obtained on all grids finer than h r = A s /2, h z = A z /10 (i.e., finer 
than the coarsest of the previous grids), and all domains larger or equal than 3 max = r max = 20, are almost 
identical. We do not plot these solutions as they are very close to one another, we rather summarize the 
results of our computations in Table 7.3, in which the two key quantities for each case are presented: The 
maximum value of self-focusing, defined as max z |£1(0, ^)| (i.e., the peak on the curve similar to those shown 
on Figure 7.3), and the maximum backscattering at z = 0, defined as max r \E(r, 0) — EPj C (r)| (i.e., the peak 
on the curve similar to those shown on Figure 7.4). 


Table 7.3 

Grid refinement and domain enlargement study for e = 0.04. 


•^m ax 

T max / 2max 

h z 

h r 

max. self-focusing 

max. backscattering 

20 

l 

A*/10 

X z /A 

1.0136 

0.013 

20 

l 

Az/10 

X z /8 

1.0129 

0.0128 

20 

2 

Az/10 

A z /4 

1.0135 

0.0128 

40 

1 

A,/10 

A,/4 

1.0132 

0.0127 

20 

1 

X z /20 

A z /4 

1.0124 

0.0112 

20 

1 

Xz/20 

X z /8 

1.0119 

0.0111 


From Table 7.3 we see that all values of maximum self-focusing that we have computed on different 
grids and different domains differ from one another by at most 0.17%. This indicates that for those ranges of 
parameters (grid sizes and domain sizes) that we have used the numerical solution is already “well converged.” 
The level of backscattering in all our simulations is between 1.1%- 1.3% of the incoming power, which again 
constitutes an error of only about 0.2% (relative to the maximum of the solution). One should probably 
regard the computational variant presented in the last row of Table 7.3 as the most accurate one because it 
was computed on the finest grid. The corresponding backscattering profile (for h z = A z /20, h r = X z /8) is 
shown in Figure 7.4(b). We again see that this profile is practically the same as the one from Figure 7.4(a), 
which corresponds to the grid twice as coarse in each direction. 

We now look at the convergence histories for our numerical solutions. Let us recall that the iteration 
scheme that we employ is nested. On the inner loop we solve a variable-coefficient linear equation, whereas 
on the outer loop we iterate with respect to the nonlinearity. Currently, we update the coefficient k 2 = 
fcp (l + e|E| 4 ) , i.e., make one nonlinear iteration, every ten linear iterations [i.e., in the notations of Section 4, 
M(n) = 10 in (4.4)]. In Figure 7.5 we show the convergence histories for the two cases that we have discussed 
before — those that correspond to the first and last rows of Table 7.3 (Figure 7.5(a) and Figure 7.5(b), 
respectively). 
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100 200 300 400 50 100 150 200 

iteration number (n) iteration number (n) 


(a) h z = \ z /10, h r = Az/4 (b) h z = Az/20, h r = A z /8 

Fig. 7.5. Convergence of iterations for e = 0.04, z m ax = 20, r ma , x /z m ax = 1. 

The actual quantity shown in Figures 7.5 is the maximum absolute difference between the two consecutive 
iterations. The sawtooth character of both curves is accounted for by the nested structure of the iterative 
procedure. The fast-scale decay followed by a jump back up is the convergence of linear iterations on the 
inner loop with subsequent update of k 2 . The slow-scale decay all the way up to machine zero corresponds 
to the convergence of nonlinear iterations on the outer loop. 

Figures 7.5 demonstrate the convergence 
of iterations. Besides, we notice that on a 
finer grid, see Figure 7.5(b), this convergence 
is faster (about twice as fast) than on the 
coarser one, see Figure 7.5(a). In fact, we have 
observed in different simulations that the ge- 
ometry in the r direction influences the rate 
of convergence most noticeably. The larger 
the domain size r max and/or the finer the grid 
size h r , the faster the iterations converge. As 
of yet, we do not have a rigorous explanation 
of this computational phenomenon. We can 
only assume that both refining the grid in the 
0 20 r direction and putting the boundary r = r max 

z further away somehow reduce the adverse in- 

FtC. 7.6. !-E compu ted(0, z)\ for the two-way ABC (solid) and for fluence q£ tMs re fl e cting boundary on the SO- 
the Dirichlet BC (dots). 

lution. 

As stated at the beginning of this section, a major goal of the nonlinear simulations is to compare the 
performance of the new two-way ABC against that of the traditional Dirichlet boundary condition at z = 0 
(7.7). In Figure 7.6 we compare the actual computed solutions with the two boundary conditions for the 
case that we have analyzed before: e = 0.04, ;z max = 20, r max /z max = 1, h z = A z /10, h r = A a /4. We see a 
noticeable discrepancy between the two curves. The dotted line that corresponds to the Dirichlet boundary 



33 






LU 

I 

c 

LU 


X 

to 

E 



Fig. 7.7. Same as Figure 7.5(a), with the Dirichlet boundary 
condition at z = 0. 


conditions is above the solid one, which corresponds to the two-way ABC. The extent of the aforementioned 
discrepancy is roughly equal to the level of backscattering that we have recovered previously, which is clearly 
a natural result to observe. 

We also compare the rates of convergence 
of the iterative algorithm for the two types 
of boundary conditions that we set at z = 0. 
The convergence history for the two-way ABC 
is shown in Figure 7.5(a), the convergence his- 
tory for the Dirichlet boundary conditions is 
shown in Figure 7.7. We see that the con- 
vergence with the two-way ABCs is about 
1.5 times faster that that with the Dirich- 
let boundary conditions, which presents an- 
other advantage of using the new methodol- 
ogy. Let us mention that the phenomenon 
of convergence speedup for iterative solvers 
caused by the application of highly-accurate 
nonlocal ABCs (similar to those developed in 
this paper) has been noticed previously by sev- 
eral authors, although in completely different settings primarily associated with the fluid flow computations, 
see [28]. 

We now consider the case e = 0.06, for 
which the input beam power is 90% of the 
critical power. Basically, the results have 
the same qualitative features as for the case 
e = 0.04. In particular, the convergence of it- 
erations is faster for finer grids and larger com- 
putational domains, as well as for the two-way 
ABC compared with the traditional Dirichlet 
boundary condition at 2 = 0. Moreover, we 
note that for e = 0.06 some cases with the 
Dirichlet boundary condition did not converge 
at all. 

In Figure 7.8A, we plot the on-axis am- 
plitude raised to the power 4 for the domain 
of the same size as corresponds to Figure 7.3 
(but with a finer grid). We plot this particu- 
lar quantity because on one hand, it is the one 
that controls the relative magnitude of nonlinearity, which is crucial for our study, and on the other hand it 
allows to see most clearly that the solution for 3 max = 20 has small oscillations throughout the domain, which 
are reminiscent of those seen in Figure 7.2. In order to verify that these oscillations are indeed due to the 
right boundary z = z max being placed too close, we re-ran the same simulation but with the right boundary 
located at twice the previous distance, i.e., ;z max = 40. The corresponding profile of | ^computed (0, z)\ A is 




Fig. 7.8. |£ CO mputed(0, z) I 4 for e = 0.06, h z = A z /20, h r = A z /8, 
Imax/ -rnax = 1- A -max = 20, B -max — 40. 
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shown in Figure 7.8B. but only for the half of the new range: 0 < z < 20, to make the scale the same as 
that on Figure 7.8A. From Figure 7.8B we see that in the case 2 max = 40 the little wiggles have almost 
disappeared, suggesting that this is indeed a numerical artifact, rather than a true physical phenomenon. 
Apart from the little wiggles, the two solutions seem to be identical as Figure 7.9 indicates. 

The explanation for the appearance of the 
small wiggles throughout the domain when the 
right boundary is too close is the same as in 
the linear case, namely, that in order for the 
ABC at 2 max to perform well, e|i?| 4 should be 
sufficiently small there so that k 2 rs with 
sufficient accuracy. Therefore, at higher e, one 
needs more decay in |f?|' for this approxima- 
tion to hold. On top of that, at higher pow- 
ers self-focusing is stronger, implying that \E\ 4 
would decay slower in 2 . This, in turn, means 
that we may need to take larger and larger do- 
mains at higher powers, otherwise, the quality 
of the computed solution will deteriorate. Be- 
sides, the convergence rate of our iterations 
may also be affected by the location of the 
boundary z = z m ax . For higher powers on 
those domains that we have considered it becomes prohibitively slow (if there is convergence at all). This is 
the reason why, at present, we could not go above e = 0.06. We should note, however, that besides enlarging 
the domain, changing the iterative algorithm itself to a more efficient one may alleviate the aforementioned 
problem. This issue will be studied in the future. 

The results of the grid convergence study for e = 0.06 are summarized in Table 7.4. Comparison of 
Table 7.3 with Table 7.4 shows that as the input power increases (relative to the critical power), more energy 
gets back-scattered and the self-focusing peak becomes higher, which is expected from physical considerations. 


Table 7.4 

Same as Table 7.3 with e = 0.06. 


^max 

T max / ^max 

h z 

hf> 

max. self-focusing 

max. backscattering 

20 

1 

A./10 

Aj/4 

1.0567 

0.0188 

20 

1 

As/20 

As/8 

1.0528 

0.0188 

20 

1 

As/20 

A z /16 

1.0526 

0.0188 

20 

2 

As/20 

As/8 

1.0527 

0.0188 

20 

1 

As/40 

As/8 

1.0518 

0.0179 

40 

1 

As/20 

As/8 

1.0512 

0.0173 



Fig. 7.9. |F?c Om p U te<i(0> ^)l ^ — 0.06, hz — As/20, hr — A 2 /8, 

rmax/zmax = 1. Solid line — z m ax = 20, dotted line — z ma x = 40. 


8. Discussion. In this section we briefly describe the approaches that have been used previously in the 
literature for solving similar problems, discuss the motivation behind making some particular choices when 
constructing our algorithm, present the conclusions, and outline directions for future research. 
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8.1. Previous approaches for solving the NLH. Feit and Fleck solved the NLH by splitting the 
wave into its forward- and backward-components, and solving only for the forward propagating part. Under 
this approach it was assumed that the “transverse variation in [k] is sufficiently small.” As for backscattering, 
their algorithm “removes power that cannot propagate in the forward direction without accounting explicitly 
to where it goes” [8]. Akhmediev and collaborators [1,2] solved an initial- value problem which can be viewed 
as a “modified” NLH. However, they neglected the ip zz term, as well as backscattering. 

In contrast to the aforementioned approaches, in this paper we solve the Helmholtz equation as a true 
“unabridged” boundary value problem. By doing that, we can account correctly for the backscattering, 
without introducing any ad-hoc assumptions, the validity of which is unclear. 

8.2. Discontinuity at the interface z = 0. In the current study we consider the simplest possible 
model for the interface z = 0, where we assume that this interface is non-reflecting, i.e., the wavenumber k is 
continuous across z = 0 (Section 2.1). From the standpoint of physics this is, of course, not necessarily true. 
For example, an incoming laser beam traveling through air which impinges on a water interface would be 
partially reflected, due to the difference in the (linear) index of refraction between air and water. The easiest 
way to incorporate the discontinuity in k at z = 0 into the model would be to do that already for the linear 
constant-coefficient equation (4.4) in the framework of the iteration scheme, as we do all other boundary 
conditions. After the transverse Fourier transform, we obtain a collection of one-dimensional Helmholtz 
equations. For each of the latter, the application of the standard elliptic interface conditions, which for the 
second-order equations are the continuity of the solution and its flux across the interface, yields the standard 
expressions for the reflection and transmission coefficients, once the incoming wave is given. If we want to 
use the transmitted wave (i.e., already past the interface) as the primary data for the problem, the same 
expressions will yield the amount of reflections and the original incoming signal. Moreover, they will also 
apply to treating the possible reflection of the backscattered waves by the interface z = 0. 

8.3. Nonlinear iterations. The primary motivation behind our choice of the nonlinear iteration 
scheme (see Section 4) was its simplicity. We note that equations (4.1), (4.2) have been obtained by simply 
freezing the nonlinear term rather than differentiating it in the sense of Frechet. For complex- valued solu- 
tions E (which is the case in our study) the nonlinearity in equation (3.2) is obviously non-differentiable 
and consequently, the direct implementation of the Newton’s method is not possible. As, however, been 
mentioned by Bayliss [3], Newton-type iterations may still apply to equation (3.2) if it is solved separately 
for the real and imaginary components of E. We did not try to implement this idea in the current study. 
We acknowledge, however, that among the different parts of our algorithm the nonlinear iteration scheme 
is apparently the primary candidate for improvements in order to achieve convergence with higher input 
power, i.e., for larger e. 

8.4. Linear solver. The solver that we employ for the variable-coefficient linear Helmholtz equation is 
also iterative and fits as the inner loop of the overall nonlinear solver. This choice is, of course, by no means 
unique. In general, one can solve the linear Helmholtz equation with variable coefficients using a variety of 
other methods, such as the Ricatti method [16]. A recent review of different approaches for solving the linear 
Helmholtz equation by Turkel can be found in [29]. We note, however, that combining a Helmholtz solver 
with global ABCs, and in particular, a two-way ABC of the type constructed in this paper, presents a rather 
difficult task, since the speed of propagation of plane waves in the z direction depends on their transverse 
wavenumber. Indeed, most of the solvers available in the literature deal with simpler boundary conditions, 
such as those of the Dirichlet type. The solver that we have constructed involves a direct inversion of the 
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constant-coefficient operator on every iteration using the separation of variables. This approach, as has been 
mentioned, is most natural for incorporating global ABCs into the model. 

8.5. Fourth-order scheme. In this study we chose a fourth-order method, rather than a conventional 
second-order one, for our simulations. The motivation behind this choice is, in fact, standard, and relies 
primarily on the possibility of having less points per wavelength and accordingly reducing the required 
overall grid dimension for a given level of accuracy. Besides, our numerical simulations corroborate that the 
extent of backscattering in the model that we study is indeed small. In the cases like that, i.e., when the 
interesting phenomenon is small in magnitude compared to the background, it is generally acknowledged 
that higher-order methods perform better than lower-order ones. 

We note in this connection that the construction of one-way discrete Helmholtz equations and radiation 
ABCs for a second-order scheme would be conceptually the same as the one described in Section 6 but 
substantially less cumbersome in both derivation and implementation, as it would not require taking care 
of an extra pair of evanescent waves. However, having a higher order method justifies, in our opinion, the 
additional work invested in obtaining the more sophisticated ABCs. 

8.6. Conclusions. Summarizing, we say that in the current paper we have developed and implemented 
a fourth-order finite-difference method for solving the nonlinear scalar Helmholtz equation that accounts for 
the phenomena of self-focusing and backscattering. The method is supplemented by the highly-accurate 
global ABCs that make the external artificial boundaries fully transparent for all outgoing waves (including 
the backscattered waves) and at the same time are capable of correctly prescribing the incoming signal at 
the outer boundary of the computational domain. To the best of our knowledge this is the first attempt ever 
of constructing global ABCs that possess the foregoing two-way capability. 

The fourth-order grid convergence of the method has been directly verified by solving model linear 
problems. In the presence of backscattering, the new method clearly outperforms a traditional technique 
based on the Dirichlet boundary condition. We have also conducted a comprehensive experimental study of 
the nonlinear case in the regime where the input power is below the critical one for blowup. Similarly to the 
linear case, this study corroborates the convergence of the method and its superiority over the traditional 
approach. 

The new method allows for a systematic quantitative study of backscattering in nonlinear self-focusing. 
To the best of our knowledge, this is the first study that allows, for example, to calculate the actual extent 
of backscattering, its dependence on the input power, etc. As has been mentioned, the new extended 
capabilities are accounted for by the fact that, unlike previous studies, we solve the NLH as a true nonlinear 
boundary value problem, without introducing any simplifying assumptions on the continuous level prior to 
the discretization. Therefore, the only error that we are actually left with is the truncation error associated 
with the discrete approximation of derivatives. 

8.7. Future work. In this paper we have developed a. new numerical methodology for solving the true 
boundary value problem for the NLH. We believe that our approach can be extended to address various 
other issues that are not covered by the present study. For example, it is interesting to conduct a systematic 
comparison of NLH simulations with the corresponding NLS simulations. Such a comparison would enhance 
our understanding on the role of nonparaxiality and backscattering. It is also interesting to compare our 
NLH simulations with the earlier approaches for solving the NLH, which did not treat the NLH as a true 
boundary-value problem. In addition, future studies should attempt to go above the critical power for 
blowup. If successful, this would provide a strong support for the current belief that there is no blowup in 
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the presence of nonparaxiality. 

In this study we have primarily focused on the NLH which corresponds to the critical NLS. However, our 
numerical approach can be applied for both subcritical NLS (e.g., calculating the amount of backscattering 
for solitons), as well as the supercritical case. 

We finally note that the nonlocal homogeneous radiation ABC at z = z max , as well as the nonlocal non- 
homogeneous two-way ABC at z = 0, can be cast into the general framework of pseudo-differential boundary 
equations and projection operators of Calderon’s type (the Calderon equation in the case of the two-way 
ABC will be non-homogeneous as well) and the difference potentials method by Ryaben’kii, see [5,18,21-24]. 
This, in particular, may allow considering curvilinear outer boundaries if necessary, as opposed to only 
linear boundaries considered in the current study. Besides, such a reformulation will be generally useful from 
the standpoint of understanding the fundamental connections between global ABCs of different types that 
appear in the scientific computing literature. 
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